Health Indicators Construction for Damage Level Assessment in Bearing Diagnostics: A Proposal of an Energetic Approach Based on Envelope Analysis

Predictive maintenance strategies are established in the industrial context on account of their benefits in terms of costs abatement and machine failures reduction. Among the available techniques, vibration-based condition monitoring (VBCM) has notably been applied in many bearing fault detection problems. The health indicators construction is a central issue for VBCM, since these features provide the necessary information to assess the current machine condition. However, the relation between vibration data and its sources intimately related to bearing damage is not effortlessly definable from a diagnostic perspective. This study discloses a diagnostic investigation performed both on the vibration signal and on the contact pressure signal that is supposed to be one of main forcing terms in the dynamic equilibrium of the damaged bearing. Envelope analysis and spectral kurtosis (SK) are applied to extract and compare diagnostic features from both signals, referring to the Case Western Reserve University (CWRU) case-study. Namely, health indicators are constructed by means of physical considerations based on the effect of faults on the signal power contents. These indicators show to be promising not only for damage detection but, also, for damage severity assessment. Moreover, they provide an invaluable reading key of the link occurring between the contact pressure path and the vibration response.


Introduction
Many studies in recent years have been devoted to the development of diagnostics techniques able to assess machinery health conditions from vibration data. The past twenty years have seen increasingly rapid advances in this field, hence the needs of the growing interest in condition-based monitoring (CBM) have met an ever wider range of monitoring strategies [1]. CBM plays a key role in the industrial context, where machines run for long periods and the long-term economic advantages related to predictive maintenance are significant, especially if compared to preventive maintenance ones [2]. Indeed, even though predictive maintenance requires additional instrumentation, such as sensors, acquisition systems and data elaboration tools, machine conditions are known at any time by means of health indicators properly constructed from the available data. Therefore, the remaining useful life (RUL) can be estimated on the basis of the current conditions, thus saving costs and time with respect to time-based preventive maintenance strategies that mostly rely on a priori time-to-failure assessments. In these cases, some components such as rolling element bearings may suffer a large the vibration data in order to evaluate a correlation between the diagnostic content of the two signals. Given such considerations, the damage detection analysis is conducted on vibration data coming from the Case Western Reserve University (CWRU) Bearing Data Center [55], which has become a benchmark dataset for bearing diagnostics [12,38,[56][57][58][59][60].
In this paper, it is proven that the energetic health indicators not only give insight into damage detectability and damage severity estimation, but they also relate the diagnostic contents of two different signals. Further, this seems to corroborate the hypothesis that the two signals are, to some extent, connected and that the energetic indicator is an evaluable reading key of such a link. Nevertheless, this preliminary study needs to be validated by a more refined statistical analysis considering a significant amount of data.

Materials and Methods
In this work, the CWRU dataset was investigated to gain an understanding of the interaction between contact pressure and vibrations in bearings with localized faults. Diagnostic parameters were specifically constructed to rate damage severity and to enable a comparison between the dynamic response and one of its main sources of excitation. In this regard, well-established methodologies such as SK-based band-pass filtering and envelope analysis were employed. Namely, the first contributed to the enhancement of the nonstationary part, whereas the second took action in the signal demodulation. The Hilbert transform adopted for signal demodulation is reported in Equation (1): where x(τ) is the time signal, and x(t) is the Hilbert transform of the signal.
Signal processing was carried out using a MATLAB ® script specifically developed for this purpose. All the details related to the signal-processing parameters will be discussed in the following sections. The pressure signal was simulated for the bearing under consideration by means of a non-Hertzian contact model implemented in a MATLAB ® code. Namely, the code was built by following the main assumptions reported in the work of Marmo et al. [54] for the numerical solution of contact problems. The Ball Passing Frequency on the Outer race (BPFO), on the Inner race (BPFI) and the fundamental train frequency (FTF) of Equations (2)-(4) can be expressed as functions of the rotation frequency f r , the number of rolling elements n,the diameter of the rolling elements d and the load angle φ.

Case Western Reserve University Bearing Data Set (CWRU)
CWRU test apparatus (Figure 1) consists of a 2-hp electric motor, a torque transducer/encoder, dynamometer and control electronics. The motor shaft is supported by test ball bearings, indicated as fan end (FE) and drive end (DE) bearings. The specifications of the tested components are reported in Table 1.  Localized faults were introduced by means of an electro-discharge machining (EDM) on the inner race, outer race and balls with fault diameters varying from 0.007 to 0.040 in (0.18-1.02 mm). Defect depth was 0.011 in (0.28 mm) for diameters between 0.007 and 0.021 in, whereas it was 0.050 in (1.27 mm) for diameters of 0.028 and 0.040 in. The details of all the fault specifications are reported in [55]. In the present work, the case of an inner race damage introduced in the DE bearing was analyzed. In these circumstances, the damages with diameters of 0.007, 0.014 and 0.021 in (0.18, 0.35 and 0.53 mm) were considered, since they referred to the same bearing type (SKF 6205-2RS JEM). Tests were carried out for each damage level with the absorbed powers of 0, 1, 2 and 3 hp (0, 0.735, 1.470 and 2.205 kW). In these tests, the running speed N ranged between 1721 and 1796 rpm. Even though SI units (International System of Units) were preferred by the author of this work, they also chose to preserve coherence with the original work of the CWRU and with the extensive literature that relates to them.
Vibration data were collected using accelerometers fixed to the motor housing by means of a magnetic base at the 12 o'clock position (opposite to the load zone-6 o'clock position). The radial load, as pointed out by Smith and Randall [38], is mostly given by the static gravitational load. Baseline and fault-bearing signals were acquired using a 16-channel DAT recorder. Digital data handled in this study were sampled at 48 kHz (Table 2). Samples were preprocessed by extracting parts with equal time durations in order to take advantage of the automatable procedure adopted in this study. In this context, it was verified that the frequency resolution provided by the choice of the number of samples M was able to correctly represent the fast Fourier-transform (FFT) spectra. Furthermore, a reasonable number of fault impulses (more than 200) was contained in the time duration T. Localized faults were introduced by means of an electro-discharge machining (EDM) on the inner race, outer race and balls with fault diameters varying from 0.007 to 0.040 in (0.18-1.02 mm). Defect depth was 0.011 in (0.28 mm) for diameters between 0.007 and 0.021 in, whereas it was 0.050 in (1.27 mm) for diameters of 0.028 and 0.040 in. The details of all the fault specifications are reported in [55]. In the present work, the case of an inner race damage introduced in the DE bearing was analyzed. In these circumstances, the damages with diameters of 0.007, 0.014 and 0.021 in (0.18, 0.35 and 0.53 mm) were considered, since they referred to the same bearing type (SKF 6205-2RS JEM). Tests were carried out for each damage level with the absorbed powers of 0, 1, 2 and 3 hp (0, 0.735, 1.470 and 2.205 kW). In these tests, the running speed N ranged between 1721 and 1796 rpm. Even though SI units (International System of Units) were preferred by the author of this work, they also chose to preserve coherence with the original work of the CWRU and with the extensive literature that relates to them.
Vibration data were collected using accelerometers fixed to the motor housing by means of a magnetic base at the 12 o'clock position (opposite to the load zone-6 o'clock position). The radial load, as pointed out by Smith and Randall [38], is mostly given by the static gravitational load. Baseline and fault-bearing signals were acquired using a 16-channel DAT recorder. Digital data handled in this study were sampled at 48 kHz (Table 2). Samples were preprocessed by extracting parts with equal time durations in order to take advantage of the automatable procedure adopted in this study. In this context, it was verified that the frequency resolution provided by the choice of the number of samples M was able to correctly represent the fast Fourier-transform (FFT) spectra. Furthermore, a reasonable number of fault impulses (more than 200) was contained in the time duration T.

Non-Hertzian Contact Model
In order to obtain accurate pressure distributions, even in presence of a defect, a frictionless non-Hertzian numerical 3D surface contact model was used. This formulation was adopted since it is able to describe pressure distributions resulting from nonconstant curvatures and sudden geometry variations. Then, the discontinuities represented by defects are analyzed by means of this model. The contact conditions can be expressed in the so-called Hertz-Signorini-Moreau problem [61][62][63].
The first condition enforces that no interpenetration can occur between the contact bodies, and therefore, the gap function h, which measures the distances between the surfaces, can only be positive or equal to 0 in the contact area. In Equation (5), and in the following, bold lowercase variables represent vectors, while bold uppercase variables represent matrices. The second condition imposes that the contact is nonadhesive, and therefore, no tension force can be present in the contact area formulated from the normal stress: where t is the traction force vector, n is the normal direction to the surface and p n = −σ n . The third condition enforces that the normal pressures can only be different from 0 inside the contact area where h = 0 and null everywhere else. The gap function h is expressed as where h 0 is the indentation between the profiles imposed as a rigid body motion, g is the initial separation of the contacting surfaces and represents its topography, while δ represents the elastic deformation of the surfaces due to the applied normal pressure p n and can be expressed as [64] δ = C · p n (8) where C is a matrix of the influence coefficients, which introduces the elasticity of the contacting surfaces. Its components C i,j (i, j = 0, 1, . . . , N) relate the displacement δ i at a point i due to the application of a unit pressure at point j.
The pressure distribution is estimated on the triangular grid-forming linearly varying pyramidal pressure elements, the definition of the influence coefficients C i,j is detailed in [65,66] and the solution procedure, which needs two iterative loops-one to eliminate tension forces and one to obtain the correct initial indentation h 0 -is described in [54]. The computational model to which this contact algorithm will be applied is visible in Figure 2, alongside a detail of the representation of the geometry of the defect that will be introduced on the inner race for the reference size of 0.014 in (0.35 mm).

Non-Hertzian Contact Model
In order to obtain accurate pressure distributions, even in presence of a defect, a frictionless non-Hertzian numerical 3D surface contact model was used. This formulation was adopted since it is able to describe pressure distributions resulting from nonconstant curvatures and sudden geometry variations. Then, the discontinuities represented by defects are analyzed by means of this model. The contact conditions can be expressed in the so-called Hertz-Signorini-Moreau problem [61][62][63].
The first condition enforces that no interpenetration can occur between the contact bodies, and therefore, the gap function , which measures the distances between the surfaces, can only be positive or equal to in the contact area. In Equation (5), and in the following, bold lowercase variables represent vectors, while bold uppercase variables represent matrices. The second condition imposes that the contact is nonadhesive, and therefore, no tension force can be present in the contact area formulated from the normal stress: where is the traction force vector, is the normal direction to the surface and = − .
The third condition enforces that the normal pressures can only be different from inside the contact area where = and null everywhere else. The gap function is expressed as where ℎ is the indentation between the profiles imposed as a rigid body motion, is the initial separation of the contacting surfaces and represents its topography, while represents the elastic deformation of the surfaces due to the applied normal pressure and can be expressed as [64] = ⋅ where is a matrix of the influence coefficients, which introduces the elasticity of the contacting surfaces. Its components , ( , = 0,1, … , ) relate the displacement at a point due to the application of a unit pressure at point .
The pressure distribution is estimated on the triangular grid-forming linearly varying pyramidal pressure elements, the definition of the influence coefficients , is detailed in [65,66] and the solution procedure, which needs two iterative loops-one to eliminate tension forces and one to obtain the correct initial indentation ℎ -is described in [54]. The computational model to which this contact algorithm will be applied is visible in Figure 2, alongside a detail of the representation of the geometry of the defect that will be introduced on the inner race for the reference size of 0.014 in (0.35 mm).  Although only gravitational loads were supposed, since no direct information regarding the actual load applied to the damaged bearing was given and the details regarding the apparatus, the motor and the related weights, three load levels were hypothesized to have more data to extrapolate from. Those load levels were inferred from the bearing specifics and, indeed, from the manufacturer of the bearing; the fatigue load limit is given as P u = 335 N, and therefore, the force F to be applied to the contact model was chosen as F 1 = 1 8 P u , F 2 = 1 4 P u and F 3 = 1 2 P u . The experimental activity was conducted in such a way that only the gravitational load acted on the bearings, but the information about weights and motor specifications were not available. Therefore, the numerical signal was simulated by hypothesizing three load levels below the fatigue threshold and quite different among them. The first condition encompassed the assumption of an infinite life design, whereas the second enabled discussions on the variability of the diagnostic content with respect to the load level, as well as to damage the severity. Examples of the obtained pressure distributions for those force levels and different defect sizes are visible in Figure 3. Indeed, the pattern of the pressure distribution is Hertzian only when no defect is present in the contact area, while introducing the defect in the middle of the contact area introduces further overpressures, due to the sharp edges of the defect and the lack of material underneath the contact area. Furthermore, since the dimension of the damage is larger than the contact area, it produces two separate contact patterns in which the maximum pressure values are around three times higher than under normal conditions. The shape and the curvatures of the contact surfaces are simulated in the MATLAB ® code by referring to the bearing geometry specifications. The shapes and the sizes of the damages, for which an example is reported in Figure 2b, retrace the experimental data of Table 2. This model was then being used to extract the maximum contact pressure as a function of the rotation angle, numerically simulating around 100 revolutions in order to have a statistically significant sample, since the BPFI is not an integer with respect to a single revolution. Computations were carried out by simulating the relative motion between the race and the rolling element. No defect evolution due to plastic deformation was taken into account, since low loads and low numbers of cycles were simulated. The obtained time histories for three different load levels and damage dimensions are shown as examples in Figure 4 for the first few revolutions of the inner race.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 25 Although only gravitational loads were supposed, since no direct information regarding the actual load applied to the damaged bearing was given and the details regarding the apparatus, the motor and the related weights, three load levels were hypothesized to have more data to extrapolate from. Those load levels were inferred from the bearing specifics and, indeed, from the manufacturer of the bearing; the fatigue load limit is given as = 335 N, and therefore, the force to be applied to the contact model was chosen as = , = and = . The experimental activity was conducted in such a way that only the gravitational load acted on the bearings, but the information about weights and motor specifications were not available. Therefore, the numerical signal was simulated by hypothesizing three load levels below the fatigue threshold and quite different among them. The first condition encompassed the assumption of an infinite life design, whereas the second enabled discussions on the variability of the diagnostic content with respect to the load level, as well as to damage the severity. Examples of the obtained pressure distributions for those force levels and different defect sizes are visible in Figure 3. Indeed, the pattern of the pressure distribution is Hertzian only when no defect is present in the contact area, while introducing the defect in the middle of the contact area introduces further overpressures, due to the sharp edges of the defect and the lack of material underneath the contact area. Furthermore, since the dimension of the damage is larger than the contact area, it produces two separate contact patterns in which the maximum pressure values are around three times higher than under normal conditions. The shape and the curvatures of the contact surfaces are simulated in the MATLAB ® code by referring to the bearing geometry specifications. The shapes and the sizes of the damages, for which an example is reported in Figure  2b, retrace the experimental data of Table 2. This model was then being used to extract the maximum contact pressure as a function of the rotation angle, numerically simulating around 100 revolutions in order to have a statistically significant sample, since the BPFI is not an integer with respect to a single revolution. Computations were carried out by simulating the relative motion between the race and the rolling element. No defect evolution due to plastic deformation was taken into account, since low loads and low numbers of cycles were simulated. The obtained time histories for three different load levels and damage dimensions are shown as examples in Figure 4 for the first few revolutions of the inner race.    Table 3 reports the main features of the simulated contact pressure signals. The initial sampling frequency was rather high (74 kHz), given the fact that it was obtained numerically. Hence, the pressure signal was numerically resampled by means of interpolation and an antialiasing filter in order to make the damage detection analysis comparable with the experimental investigation. To the same end, numerical samples were cut by using M points.

Damage Detection Analyisis
The damage detection analysis was carried out by taking advantage of the well-established concepts presented in the introduction to this work. Namely, different diagnostics features were extracted from both signals. Although the greatest part of these features enabled damage detectability, the fault severity assessment was not straightforward. Furthermore, it was a hardly identifiable correlation between the two signals in a diagnostic sense. For this reason, the authors proposed two health indicators to perform these tasks. Indeed, the dedicated features were constructed in an attempt to investigate direct correlations. Additionally, the conjectured sensitivity to the damage level was endorsed by the experimental validation. Vibration data coming from FE healthy bearing were analyzed as well. Prior to introducing the proposed indicators, the standard features and the methodology for their extraction are presented in the following.   Table 3 reports the main features of the simulated contact pressure signals. The initial sampling frequency was rather high (74 kHz), given the fact that it was obtained numerically. Hence, the pressure signal was numerically resampled by means of interpolation and an antialiasing filter in order to make the damage detection analysis comparable with the experimental investigation. To the same end, numerical samples were cut by using M points.

Damage Detection Analyisis
The damage detection analysis was carried out by taking advantage of the well-established concepts presented in the introduction to this work. Namely, different diagnostics features were extracted from both signals. Although the greatest part of these features enabled damage detectability, the fault severity assessment was not straightforward. Furthermore, it was a hardly identifiable correlation between the two signals in a diagnostic sense. For this reason, the authors proposed two health indicators to perform these tasks. Indeed, the dedicated features were constructed in an attempt to investigate direct correlations. Additionally, the conjectured sensitivity to the damage level was endorsed by the experimental validation. Vibration data coming from FE healthy bearing were analyzed as well. Prior to introducing the proposed indicators, the standard features and the methodology for their extraction are presented in the following. Time domain and frequency domain features were investigated considering their variability linked to the damage severity and absorbed power (for the vibration signal) or applied load (for the pressure signal). Root mean square value (RMS), crest factor (CF) and kurtosis (K) are reported in Equations (9)- (11).
Appl. Sci. 2020, 10, 8131 where T is the period, x(t) is the time signal, x peak is the peak value of the signal, E() is the expected value operator, µ is the mean value and σ is the standard deviation. Kurtosis thus computed is expected to take the value of three for gaussian signals. Frequency domain features were extracted from FFT raw spectra ( Figures 5-7). In particular, the mean contribution, the rotating contribution and the fault contribution, respectively corresponding to the 0-Hz harmonic, 1xN harmonic and 1xBPFI harmonic were individually analyzed.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 25 where is the period, ( ) is the time signal, is the peak value of the signal, () is the expected value operator, is the mean value and is the standard deviation. Kurtosis thus computed is expected to take the value of three for gaussian signals.
Frequency domain features were extracted from FFT raw spectra ( Figures 5-7). In particular, the mean contribution, the rotating contribution and the fault contribution, respectively corresponding to the 0-Hz harmonic, 1xN harmonic and 1xBPFI harmonic were individually analyzed.
where is the period, ( ) is the time signal, is the peak value of the signal, () is the expected value operator, is the mean value and is the standard deviation. Kurtosis thus computed is expected to take the value of three for gaussian signals.
Frequency domain features were extracted from FFT raw spectra ( Figures 5-7). In particular, the mean contribution, the rotating contribution and the fault contribution, respectively corresponding to the 0-Hz harmonic, 1xN harmonic and 1xBPFI harmonic were individually analyzed.   SK and FK were employed for the purpose of enhancing the bearing signal in the vibration data and isolating the impulsive part in the pressure signal. As mentioned earlier, these methods are particularly useful, though many alternatives are available in the literature. Indeed, the potential filtering bands can be underlined by comparing power spectral densities (PSD) of undamaged and damaged bearings [1]. Nevertheless, a major problem with this method is that the machine signature needs always to be known. Besides, wavelet denoising [67][68][69] works by taking advantage of wavelet transform and deleting the wavelet contributions that go beyond a defined threshold. However, there are certain drawbacks associated with the use of this methodology, since the result strongly depends on the choice of the mother wavelet and of the thresholding method. SK (Equation (12)) relies on the use of short-time Fourier transform (STFT) of Equation (13) to give a measure of the signal impulsiveness contained in frequency bands being examined.
where 〈•〉 is the time-averaging operator, ( , ) is the time-frequency transform of the signal ( ), ( − ) is a time window and is the imaginary unit.
The time-frequency representation of the signal is obtained by means of STFT, which, differently from FFT, is capable of localizing the harmonic contribution in time, since it is essentially a FFT that slides in time thanks to the window ( − ). SK is expected to be null for complex signals ( , ) originating from gaussian noises, whereas high values of SK suggest nonstationarity or non-gaussian properties. The size of the chosen window influences the SK computation, since STFT is affected by the well-known indetermination principle. Namely, the more accurate the time representation (short window), the worse the frequency resolution. For this reason, FK (Figures 8-10) represents SK as a SK and FK were employed for the purpose of enhancing the bearing signal in the vibration data and isolating the impulsive part in the pressure signal. As mentioned earlier, these methods are particularly useful, though many alternatives are available in the literature. Indeed, the potential filtering bands can be underlined by comparing power spectral densities (PSD) of undamaged and damaged bearings [1]. Nevertheless, a major problem with this method is that the machine signature needs always to be known. Besides, wavelet denoising [67][68][69] works by taking advantage of wavelet transform and deleting the wavelet contributions that go beyond a defined threshold. However, there are certain drawbacks associated with the use of this methodology, since the result strongly depends on the choice of the mother wavelet and of the thresholding method. SK (Equation (12)) relies on the use of short-time Fourier transform (STFT) of Equation (13) to give a measure of the signal impulsiveness contained in frequency bands being examined.
where · is the time-averaging operator, X( f , t) is the time-frequency transform of the signal x(t), w(t − τ) is a time window and j is the imaginary unit.
The time-frequency representation of the signal is obtained by means of STFT, which, differently from FFT, is capable of localizing the harmonic contribution in time, since it is essentially a FFT that slides in time thanks to the window w(t − τ). SK is expected to be null for complex signals X( f , τ) originating from gaussian noises, whereas high values of SK suggest nonstationarity or non-gaussian properties. The size of the chosen window influences the SK computation, since STFT is affected by the well-known indetermination principle. Namely, the more accurate the time representation (short window), the worse the frequency resolution. For this reason, FK (Figures 8-10) represents SK as a function of the window lengths. Therefore, FK provides information not only on the optimal demodulation band but, also, on the choice of the window length to be applied for the SK computation. Then, FK acts with a simplified algorithm, which avoids complete STFT computation.
Once signals were band pass-filtered in the frequency bands suggested by SK, the envelope analysis was applied. The modulus of the analytic signal was then investigated in the frequency domain ( Figure 11). The xBPFI harmonic contributions were particularly examined. Unlike raw spectra, xBPFI harmonics were fairly identifiable even in the FE signal ( Figure 12).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 25 function of the window lengths. Therefore, FK provides information not only on the optimal demodulation band but, also, on the choice of the window length to be applied for the SK computation. Then, FK acts with a simplified algorithm, which avoids complete STFT computation.   function of the window lengths. Therefore, FK provides information not only on the optimal demodulation band but, also, on the choice of the window length to be applied for the SK computation. Then, FK acts with a simplified algorithm, which avoids complete STFT computation.    Once signals were band pass-filtered in the frequency bands suggested by SK, the envelope analysis was applied. The modulus of the analytic signal was then investigated in the frequency domain ( Figure 11). The xBPFI harmonic contributions were particularly examined. Unlike raw spectra, xBPFI harmonics were fairly identifiable even in the FE signal ( Figure 12).  Once signals were band pass-filtered in the frequency bands suggested by SK, the envelope analysis was applied. The modulus of the analytic signal was then investigated in the frequency domain ( Figure 11). The xBPFI harmonic contributions were particularly examined. Unlike raw spectra, xBPFI harmonics were fairly identifiable even in the FE signal ( Figure 12).  In red, the first ten xBPFI harmonics.
It is important to underline that the filtering process reduced the range in which the signals fluctuated. Moreover, this effect was different for each signal, given different filtering bands. Therefore, damage detectability cannot purely be assessed on the basis of the absolute value of the harmonic contributions but, rather, on their identifiability. Envelope spectrum was expected to show xBPFI harmonics, with sidebands spaced at the shaft rotation frequency. Clearly, this was mostly visible in the pressure signal, where no disturbance sources (noise or discrete harmonic components) Figure 12. Envelope signal and envelope spectrum. Fan end (FE) signal with damage on DE bearing. In red, the first ten xBPFI harmonics.
It is important to underline that the filtering process reduced the range in which the signals fluctuated. Moreover, this effect was different for each signal, given different filtering bands. Therefore, damage detectability cannot purely be assessed on the basis of the absolute value of the harmonic contributions but, rather, on their identifiability. Envelope spectrum was expected to show xBPFI harmonics, with sidebands spaced at the shaft rotation frequency. Clearly, this was mostly visible in the pressure signal, where no disturbance sources (noise or discrete harmonic components) were introduced.
The xBPFI harmonic contributions were analyzed for both signal envelopes and for all the damage levels and loads. However, none of these features for the enveloped and nonenveloped signals enabled a direct comparison between the diagnostic contents of the pressure and of the vibration signal, since these quantities, though linked, were intrinsically different. In an attempt to overcome this limitation, the novel step was introduced.

Proposed Health Indicators
The leading hypothesis behind the following argument is that the informative content on the damage severity is hidden beneath the energy related to the fault harmonics in the enveloped signals. To this purpose, the authors introduced the normalized power spectral density (NPSD) [42] obtained from the ratio between the PSD and the mean power of the signal, which is RMS 2 (Equations (14)-(16)).
where E x is the signal energy.
One of the properties of the NPSD is that its integration in the frequency domain is unitary (Equation (17)).
Essentially, this is due to the fact that the integration of NPSD in a given frequency band results in the fraction of the signal power contained in such a band. Therefore, the health indicator (HI) proposed in Equation (18) gives an insight into the portion of power related to a frequency range with central frequency f c and bandwidth ∆ f .
Established that the conjunction between the SK and envelope analysis provided proper damage detectability into processed spectra, NPSD and HI were expected to give energetic information purely related to the fault, as long as they are evaluated in xBPFI harmonics. Indeed, thanks to NPSD and HI, the contribution given to the signal power solely imputable to the damage was investigated. This provided the remarkable possibility of isolating the contribution given to a global signal feature (mean power) by a very specific event, which is damage. Additionally, the physical significance of these features made the authors assume two interesting aspects. First, it was reasonable to hypothesize that the greater the defect size was, the more accentuated the relative energetic contribution given by the fault to the extracted signal. Hence, NPSD and HI were expected to be sensitive to damage severity. Secondly, the physical units of these quantities permitted a direct comparison between signals of different natures. Actually, the indicators referred to energetic quantities with respect to the total signal energy, which is different for each type of signal. Consequently, relative energetic features rather than absolute were adopted to compare different physical information in the pursuit of a correlation.
The described features were insofar evaluated in the xBPFI harmonics (Table 4) of the accelerometric and pressure signals, and the sensitivity to the load and damage level was explored. Eventually, the correlation coefficient r (Equation (19)) was considered as a measure of the linear dependence between the vibration and pressure features.
where P is the number of observations, FP is the feature extracted from the pressure signal, µ FP and σ FP are the mean and the standard deviation of FP, FV is the feature extracted from the vibration signal and µ FV and σ FV are the mean and the standard deviation of FV.

Methodology Application
As discussed previously, the aim of the study was to analyze the diagnostic content of two signals. The first is the vibration signal, which is commonly processed to extract information about bearings and machine health in a CBM perspective. The second is the contact pressure that arises between balls and race, since it is supposed to be one of the main sources of excitation in the dynamic equilibrium of the system. Therefore, the features presented in Section 2 are evaluated by means of the earlier described methodologies with the intention of comparing their diagnostic potential with that of the proposed health indicators. In most cases, although contact pressure analyses suggested smooth trends, vibration features showed to be hardly connectable to damage level and loads. Even more, a signals correlation was barely identifiable. Conversely, NPSD and HI showed to be appreciably sensitive to damage severity, and they paved the way for a direct comparison. In this section, the results of the application of the standard features and of the proposed indicators are presented. Figure 13 shows the trend of the time features. DE and FE data showed similar patterns, but the effect of the damage level was not always clearly definable. Conversely, peak values and RMS of the contact pressure stem directly from the applied load and from the damage level. Indeed, the mean part of the pressure signal is influenced by the load, whereas the overpressure (Figure 4) comes from the effects of the defect size.

Standard Features
The frequency features of Figure 14 showed a similar behavior. The nonenveloped pressure signal manifested a smooth trend fairly suitable for damage level detection. Instead, the vibration data did not unequivocally show this attitude. Once more, the extracted features revealed were hardly relatable to the damage level and the diagnostic contents of the two physical signals placed on different levels.
Following the concepts reported in Section 2.3, signals were band pass-filtered in the frequency bands computed by means of FK. Most of the pressure signals revealed the maximum SK in a frequency band upper-circumscribed by the Nyquist limit. According to the authors, this was unavoidable, since the pressure signal showed very narrow peaks that resulted in wide frequency domain contents. Hence, most of the impulsivity was retained in the highest representable frequency bands. Differently, in the vibration signal, SK attempted to emphasize the bands where high-frequency structural resonances were excited by nonstationary phenomena, which were impacts of the rolling elements.
presented. Figure 13 shows the trend of the time features. DE and FE data showed similar patterns, but the effect of the damage level was not always clearly definable. Conversely, peak values and RMS of the contact pressure stem directly from the applied load and from the damage level. Indeed, the mean part of the pressure signal is influenced by the load, whereas the overpressure (Figure 4) comes from the effects of the defect size. The frequency features of Figure 14 showed a similar behavior. The nonenveloped pressure signal manifested a smooth trend fairly suitable for damage level detection. Instead, the vibration data did not unequivocally show this attitude. Once more, the extracted features revealed were hardly relatable to the damage level and the diagnostic contents of the two physical signals placed on different levels. Following the concepts reported in Section 2.3, signals were band pass-filtered in the frequency bands computed by means of FK. Most of the pressure signals revealed the maximum SK in a frequency band upper-circumscribed by the Nyquist limit. According to the authors, this was unavoidable, since the pressure signal showed very narrow peaks that resulted in wide frequency domain contents. Hence, most of the impulsivity was retained in the highest representable frequency bands. Differently, in the vibration signal, SK attempted to emphasize the bands where highfrequency structural resonances were excited by nonstationary phenomena, which were impacts of the rolling elements.

Standard Features
Subsequently, the envelope was extracted for signals demodulation. As can be seen from Figure  15, low values of the harmonic contents were appreciable in the first three xBPFI contributions. However, as pointed out in Section 2, an assessment based on the absolute values of such harmonics in enveloped signals may be misleading. One of the consequences of filtering in different bands Subsequently, the envelope was extracted for signals demodulation. As can be seen from Figure 15, low values of the harmonic contents were appreciable in the first three xBPFI contributions. However, as pointed out in Section 2, an assessment based on the absolute values of such harmonics in enveloped signals may be misleading. One of the consequences of filtering in different bands according to the maximum SK criterion was to eradicate the similarity between DE and FE trends. Next, the xBPFI contents of the enveloped pressure signal were shown to be quite differently arranged with respect from the vibration analyses. The DE parameters were revealed to be not remarkably dependent on the damage level ( Figure 15). Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 25

Proposed Health Indicators: Application to the Case-Study
NPSD and HI can be interpreted by visualizing their significance on FFT spectra. Figures 16 and  17 report the example of the 0.021 in defect with 3 hp (vibration signal) and 0.021 and Pu/2 (contact pressure signal), but all the analyzed cases retraced this behavior. As long as the pressure signal is nonenveloped (Figure 16), most of the energetic content is concentrated in the low frequency range. The jumps in the cumulative NPSD point out the frequencies that mostly contribute to the signal PSD. Indeed, the very rapid increase detectable in the jumps corresponds to the harmonics in which the signal power is lumped. In particular, this occurred in the rotation frequency and xBPFI harmonics. When the pressure signal is enveloped, NPSD shows peaks corresponding to fault harmonics. This results in higher jumps in the cumulative NPSD. Raw vibration data did not exhibit this trend, since a great part of the signal power is given by resonance contributions. Conversely, it is interesting to notice that, in the enveloped vibration data (Figure 17), the energetic jumps are mostly related to the xBPFI content. This is consistent with the extensive literature supporting the strength of the envelope analysis for damage detectability [1,22,23,35,37,70]. NPSD and cumulative NPSD were computed for nonenveloped and enveloped signals considering all the loads and damage level cases.

Proposed Health Indicators: Application to the Case-Study
NPSD and HI can be interpreted by visualizing their significance on FFT spectra. Figures 16  and 17 report the example of the 0.021 in defect with 3 hp (vibration signal) and 0.021 and P u /2 (contact pressure signal), but all the analyzed cases retraced this behavior. As long as the pressure signal is nonenveloped (Figure 16), most of the energetic content is concentrated in the low frequency range. The jumps in the cumulative NPSD point out the frequencies that mostly contribute to the signal PSD. Indeed, the very rapid increase detectable in the jumps corresponds to the harmonics in which the signal power is lumped. In particular, this occurred in the rotation frequency and xBPFI harmonics. When the pressure signal is enveloped, NPSD shows peaks corresponding to fault harmonics. This results in higher jumps in the cumulative NPSD. Raw vibration data did not exhibit this trend, since a great part of the signal power is given by resonance contributions. Conversely, it is interesting to notice that, in the enveloped vibration data (Figure 17), the energetic jumps are mostly related to the xBPFI content. This is consistent with the extensive literature supporting the strength of the envelope analysis for damage detectability [1,22,23,35,37,70]. NPSD and cumulative NPSD were computed for nonenveloped and enveloped signals considering all the loads and damage level cases.
The HI proposed in Equation (18) provides a quantitative extent of the jumps, since it represents the fraction of the power contained in these gaps relative to the whole signal power. As anticipated, it was hypothesized that the higher this fraction is, the more accentuated the damage severity, since a greater part of the signal energy is relatable to the presence of the fault. It appeared that this assumption was valid for the pressure signal, in which it was found a linear trend of NPSD and HI with respect to the defect size ( Figure 18). As expected, the indicators, though showing a similar trend, exhibited lower values for the nonenveloped signal. Figure 19 shows the trend of these features for the enveloped and raw DE and FE signals. It can be observed that the damage is hardly detectable by means of these indicators without demodulation. On the other hand, in the enveloped signals, it is possible to allocate up to 20% of the signal power to the presence of localized faults. Interestingly, these values are comparable with the FE signal, since the adopted indicators are normalized with respect to the properties of each signal. Hence, in FE signals, low power fractions are related to low signal powers, resulting in values comparable with DE ones. Differently from K and CF, which have neutral physical units in the same way, NPSD and HI assume similar values in the vibration and in the pressure signal. Indeed, NPSD and HI refer to the relation between an isolated event (damage) and signal power, whereas K and CF globally measure the effect of the damage on the time domain shape.  The HI proposed in Equation (18) provides a quantitative extent of the jumps, since it represents the fraction of the power contained in these gaps relative to the whole signal power. As anticipated, it was hypothesized that the higher this fraction is, the more accentuated the damage severity, since  The HI proposed in Equation (18) provides a quantitative extent of the jumps, since it represents the fraction of the power contained in these gaps relative to the whole signal power. As anticipated, it was hypothesized that the higher this fraction is, the more accentuated the damage severity, since Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 25 a greater part of the signal energy is relatable to the presence of the fault. It appeared that this assumption was valid for the pressure signal, in which it was found a linear trend of NPSD and HI with respect to the defect size ( Figure 18). As expected, the indicators, though showing a similar trend, exhibited lower values for the nonenveloped signal.  Figure 19 shows the trend of these features for the enveloped and raw DE and FE signals. It can be observed that the damage is hardly detectable by means of these indicators without demodulation. On the other hand, in the enveloped signals, it is possible to allocate up to 20% of the signal power to the presence of localized faults. Interestingly, these values are comparable with the FE signal, since the adopted indicators are normalized with respect to the properties of each signal. Hence, in FE signals, low power fractions are related to low signal powers, resulting in values comparable with DE ones. Differently from K and CF, which have neutral physical units in the same way, NPSD and HI assume similar values in the vibration and in the pressure signal. Indeed, NPSD and HI refer to the relation between an isolated event (damage) and signal power, whereas K and CF globally measure the effect of the damage on the time domain shape.    Figure 19 shows the trend of these features for the enveloped and raw DE and FE signals. It can be observed that the damage is hardly detectable by means of these indicators without demodulation. On the other hand, in the enveloped signals, it is possible to allocate up to 20% of the signal power to the presence of localized faults. Interestingly, these values are comparable with the FE signal, since the adopted indicators are normalized with respect to the properties of each signal. Hence, in FE signals, low power fractions are related to low signal powers, resulting in values comparable with DE ones. Differently from K and CF, which have neutral physical units in the same way, NPSD and HI assume similar values in the vibration and in the pressure signal. Indeed, NPSD and HI refer to the relation between an isolated event (damage) and signal power, whereas K and CF globally measure the effect of the damage on the time domain shape.

Methodology Validation
In order to investigate the feature sensitivity to the damage severity, the indicators were analyzed by considering multiple powers or applied loads as different realizations of a given damage level. The boxplots of Figure 20 indicate the attitudes of the different features for the three damage levels: L1, L2 and L3, which correspond, respectively, to 0.007, 0.014 and 0.021 in. It appears that, for all the feature distributions, it is slightly identifiable as a clear distinction between damage levels, except for NPSD and HI, which emphasize an evident discrimination. This reinforces the original assumption taken for the construction of NPSD and HI and confirms that, actually, the power fraction merely imputable to the presence of the defect is very sensitive to damage severity. However, boxplots are employed to give a graphical representation of the data, but their significance must be consolidated by a larger amount of realizations for each damage level. except for NPSD and HI, which emphasize an evident discrimination. This reinforces the original assumption taken for the construction of NPSD and HI and confirms that, actually, the power fraction merely imputable to the presence of the defect is very sensitive to damage severity. However, boxplots are employed to give a graphical representation of the data, but their significance must be consolidated by a larger amount of realizations for each damage level. Detectability is studied by adopting the same approach even for the baseline and fault analysis performed on the vibration data, validating the capabilities of the proposed Indicators. As shown in Figures 21 and 22, all the features are able to discern faulty bearings, with the exception of the central frequency of the maximum SK band. Evidently, this feature is related to the research of nonstationarity, and no linear trend is expected between this parameter and damage severity. A clear separation between indicators for different damage levels is not identifiable for all the features. However, kurtosis (K) showed an increasing trend with well-separable distributions related to different damage levels. Detectability is studied by adopting the same approach even for the baseline and fault analysis performed on the vibration data, validating the capabilities of the proposed Indicators. As shown in Figures 21 and 22, all the features are able to discern faulty bearings, with the exception of the central frequency of the maximum SK band. Evidently, this feature is related to the research of nonstationarity, and no linear trend is expected between this parameter and damage severity. A clear separation between indicators for different damage levels is not identifiable for all the features. However, kurtosis (K) showed an increasing trend with well-separable distributions related to different damage levels.    Interestingly, the NPSD and HI were revealed to be promising features, since they fulfilled three main criteria. First, they can classify faulty bearings with respect to healthy ones. Secondly, they assume an increasing trend with respect to the damage level. Finally, their boxplots do not suffer excessive superimposition. Nevertheless, only a larger amount of data will confirm this attitude. One of the extremely attractive features of NPSD and HI is their values in the pressure signal and in the experimental signal of the damaged-bearing range in similar intervals for each fault size. Given the Interestingly, the NPSD and HI were revealed to be promising features, since they fulfilled three main criteria. First, they can classify faulty bearings with respect to healthy ones. Secondly, they assume an increasing trend with respect to the damage level. Finally, their boxplots do not suffer excessive superimposition. Nevertheless, only a larger amount of data will confirm this attitude. One of the extremely attractive features of NPSD and HI is their values in the pressure signal and in the experimental signal of the damaged-bearing range in similar intervals for each fault size. Given the fact that NPSD and HI are endowed with neutral physical units, the diagnostic correlation between contact pressure and vibration was investigated ( Figure 23).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 21 of 25 fact that NPSD and HI are endowed with neutral physical units, the diagnostic correlation between contact pressure and vibration was investigated ( Figure 23).  Figure 23 shows the results of the correlation study carried out by means of Equation (19) described in Section 3. The mean value of the damage-related power HI and of NPSD are computed for each fault level. These values constitute the and of Equation (19). The noteworthy correlation coefficients (Table 5) confirm these indicators to provide an evaluable key to link two different quantities from a diagnostic content perspective. This evidence reinforces the initial idea that one of the sources generating the dynamic response is actually related to this latter in a diagnostic sense. Moreover, it is found that, by taking into account diagnostic quantities, the possible correlation lies in the signal power. Indeed, this study highlights that the signal power strictly related to the  Figure 23 shows the results of the correlation study carried out by means of Equation (19) described in Section 3. The mean value of the damage-related power HI and of NPSD are computed for each fault level. These values constitute the FP and FV of Equation (19). The noteworthy correlation coefficients (Table 5) confirm these indicators to provide an evaluable key to link two different quantities from a diagnostic content perspective. This evidence reinforces the initial idea that one of the sources generating the dynamic response is actually related to this latter in a diagnostic sense. Moreover, it is found that, by taking into account diagnostic quantities, the possible correlation lies in the signal power. Indeed, this study highlights that the signal power strictly related to the damage is sensitive to fault dimensions.

Discussion and Conclusions
This work documents the diagnostic analysis carried out on faulty and healthy bearings of a CWRU vibration dataset. Namely, previous research contributions have drawn attention on SK for the enhancement of nonstationarity and on the envelope analysis for signal demodulation. No less research in bearing diagnostics is devoted to feature extraction methods based on data-driven models. However, in this study, the authors placed emphasis on physical phenomena underpinning the resulting experimental data, thus preferring diagnostic methodologies that do not suffer from an excessive detachment from physics-based considerations, which may help data interpretation. Moreover, these techniques are applied to evaluate the capability of damage severity assessments for different diagnostic features. To this aim, the diagnostic content was analyzed both in a numerical and in an experimental signal. The numerical signal was generated by investigating the contact pressure established between the balls and the defected race. Indeed, the contact pressure was supposed to be one of the main sources of excitation in the system dynamic response.
Therefore, envelope analysis was applied to extract features and evaluate their damage level sensitivity. The choice of the demodulation band was assisted by FK, which enabled the enhancement of the bearing signal in the vibration data, whereas, in the numerical signal, it contributed to the isolation of the impulsive phenomena.
Actually, out an inherent difference between the two signals must be pointed out. In the vibration spectrum, fault harmonics are buried due to the excitation of high-frequency structural resonances and SK acts by bringing to light the bands in which the nonstationary part is hidden. On the contrary, the numerical signal is not influenced by the structural dynamics of the whole test rig, and the peaks do not suffer from damping flattening. Therefore, FK emphasizes impulses that, however, are not completely undetectable in raw spectra.
Two health indicators were introduced. They were constructed relying on the assumption that the more severe the damage, the more the signal power strictly relatable to faults is accentuated. In the attempt of isolating this contribution, SK and envelope extraction were applied, resorting to their physical significance. The indicators showed their attitude in identifying faulty bearings with respect to healthy ones. Moreover, the results showed that they may be capable of assessing the damage severity as well. Another compelling aspect is that their physical units are neutral, paving the way for a comparison between the diagnostic contents of the source, which is the pressure signal, and the effect, represented by the vibration data. Indeed, when these signals are evaluated through the mentioned indicators, a correlation is found. Therefore, it is concluded that the energetic indicators may be proposed for damage detection and damage level assessment but, also, as reading keys of a correlation between the vibration and the pressure signal. One of the drawbacks of these indicators is that they precisely point the xBPFI harmonic contents, relying on the detectability of a significant harmonic contribution in these frequencies.
However, these analyses were carried out by considering a limited amount of data, both in terms of the fault level and applied loads. Then, future works should include more data in order to validate and give statistical significance to these considerations. Other fault types should also be investigated, together with the possible prognostic capabilities. The use of different diagnostic methodologies, such as alternatives to FK, squared envelope and data-driven models, could be investigated as well.