Ge Xin
email: xinge@bjtu.edu.cn

Jerome Antoni
email: jerome.antoni@insa-lyon.fr

Semi-automated diagnosis of bearing faults based on a hidden Markov model of the vibration signals

Keywords: Automated diagnosis, repetitive transients detection, repetitive transients extraction, bearing faults, hidden Markov model, mixture of Gaussians, cyclostationary signals, non-stationary operating conditions. Highlights

Due to its practical importance, the diagnosis of rolling element bearing has attracted constant interest in the scientific community. At the incipient stage of a failure, the measured vibration signal typically consists of a series of repetitive transients immerged in background noise. Although they are usually carried in high frequency bands due to the high stiffness of bearings, they are fairly weak compared with surrounding noise and other interfering signals. In addition, taking random slips and fluctuations into account, the transients produced by impacts are not strictly periodic but rather tend to be random cyclostationary. This makes the diagnosis of rolling element bearing quite challenging and, consequently, various signal processing techniques have been developed for either the detection, the identification or the extraction of the fault, whose combination asks for a high level of expertise of the user. The aim of this paper is to address all these objectives at once, in the same algorithm, by proposing a semi-automated method that requires the setting of only one parameter. It is rooted on a probabilistic model, in the form of a mixture of Gaussians, endowed with a hidden variable that indicates the occurrence of impacts. The method is shown to be optimal for detection in the Neyman-Pearson sense, it returns an envelope spectrum comparable to the best that can be obtained by other meanswhich often require a careful pre-filtering stepfrom which fault frequencies can be identified, and it eventually returns the fault signal from which subsequent severity/health indicators can be computed. There is almost no demand on the user's expertise (apart from setting the frequency resolution), even though the method does not address the decision part. The performance is investigated on synthetic signals and its robustness is also verified on several vibration signals measured on test-rigs. Results are found superior or at least equivalent to those of the conventional semi-automated method based on the fast kurtogram in combination with the envelope analysis.

Introduction

Rolling element bearings are common but fragile components used in various types of mechanical systems, and thus a foremost cause of machinery breakdown. Typical bearing defects are caused by crack, breakage, spall or uneven wear (pitting, scuffing, abrasion, erosion), often located on the matting surface of the inner race, the outer race or the rolling elements. As the rolling elements strike a defect, a series of vibration transients occur at a specific rate corresponding to the bearing "characteristic frequency" or "fault frequency". At the incipient stage, each transient resembles a damped impulse response with frequency content specified by the excited structural resonances. Due to the load distribution, the series of transients are additionally amplitude modulated by the passing period into and out of the load zone. Bearing fault signatures have been well-investigated and their characteristic frequencies are listed from given kinematic and geometric parameters [START_REF] Randall | Rolling element bearing diagnostics-a tutorial[END_REF][START_REF] Smith | Rolling element bearing diagnostics using the Case Western Reserve University data: A benchmark study[END_REF]. Nevertheless, bearing fault signals are often masked by numerous extraneous sources of vibration, simply referred to hereafter as "background noise", such that signal processing is necessary to recover the diagnostics information. Another difficulty arises from machines which do not operate in stationary conditions, such as wind turbines or crushers, which ask for even more advanced data processing.

The diagnosis of rolling element bearingsas understood in the present paperincludes three main steps, whose chronological order actually coincides with their degree of difficulty: 1) detection of the fault, 2) identification of the fault and 3) extraction of the fault signal. Since they correspond to different goals, they are usually addressed by different signal processing tools.

The detection of the fault is usually formulated as a problem of detecting the presence of a weak signal in background noise, ideally in the form of statistical test. In the case of bearing faults, the signal has a complex statistical structureit is actually cyclostationary, i.e. non-stationary random with periodic statistics [START_REF] Antoni | Cyclostationary modelling of rotating machine vibration signals[END_REF] which sometime makes this task challenging. This step has nourished a vast literature on a variety of methods which nowadays probably culminate to their asymptote in terms of performance. Some typical tools are time-frequency analysis [START_REF] Boashash | Time-Frequency Signal Analysis and Processing[END_REF][START_REF] Feng | Recent advances in time-frequency analysis methods for machinery fault diagnosis: a review with application examples[END_REF][START_REF] Chandra | Fault detection in rotor bearing systems using time frequency techniques[END_REF], wavelet analysis [START_REF] Ghods | Probabilistic frequency-domain discrete wavelet transform for better detection of bearing faults in induction motors[END_REF][START_REF] Shi | Intelligent bearing fault signature extraction via iterative oscillatory behavior based signal decomposition (IOBSD)[END_REF][START_REF] Yuan | Weak fault detection and health degradation monitoring using customized standard multiwavelets[END_REF], the stochastic resonance [START_REF] Lai | Weak-signal detection based on the stochastic resonance of bistableDuffing oscillator and its application in incipient fault diagnosis[END_REF][START_REF] Guo | Multi-frequency weak signal detection based on multi-segment cascaded stochastic resonance for rolling bearings[END_REF], morphological analysis [START_REF] Li | Gear fault detection using multi-scale morphological filters[END_REF], the spectral kurtosis and the fast kurtogram [START_REF] Wang | Spectral kurtosis for fault detection, diagnosis and prognostics of rotating machines: A review with applications[END_REF], the spectral correlation [START_REF] Antoni | Cyclic spectral analysis in practice[END_REF][START_REF] Antoni | Fast computation of the spectral correlation[END_REF], sparsity [START_REF] Du | Weighted low-rank sparse model via nuclear norm minimization for bearing fault detection[END_REF], etc. It must be said that in most of these methods the detection of a fault is done visually; few of them actually tackle the problem as a statistical test in terms of probabilities of detection and of false alarm.

Next, the identification of the fault mainly consists in estimating the fault frequency and associating it to a given component in the machine. The prevailing method in the modern literature is surely the squared envelope spectrum (SES). The envelope spectrum is a mean to demodulate a non-stationary signalof possible random nature, in particular cyclostationary-and identify periodic modulations related with the bearing characteristic frequencies [START_REF] Mcfadden | Vibration monitoring of rolling element bearings by the high-frequency resonance technique-a review[END_REF]. In its most frequent version, the envelope is estimated by squaring the signalor better, the modulus of the analytical signaland its Fourier spectrum is then computed. The fact that the SES contains most of the diagnostic information can be justified theoretically by modelling the bearing fault signals as cyclostationary [START_REF] Randall | The relationship between spectral correlation and envelope analysis in the diagnostics of bearing faults and other cyclostationary machine signals[END_REF]. It is noteworthy that, in the cyclostationary framework, the SES can also be introduced as a statistics for fault detection [START_REF] Antoni | Cyclic spectral analysis of rolling-element bearing signals: facts and fictions[END_REF]. However, the SES of the raw signal is rarely a good diagnostic indicator when used without pre-processing, because it is highly sensitive to the presence of background noise and other interfering components [START_REF] Randall | The relationship between spectral correlation and envelope analysis in the diagnostics of bearing faults and other cyclostationary machine signals[END_REF][START_REF] Antoni | Cyclic spectral analysis of rolling-element bearing signals: facts and fictions[END_REF][START_REF] Antoni | Differential diagnosis of gear and bearing faults[END_REF][START_REF] Abboud | Envelope analysis of rotating machine vibrations in variable speed conditions: A comprehensive treatment[END_REF][START_REF] Ming | Envelope calculation of the multi-component signal and its application to the deterministic component cancellation in bearing fault diagnosis[END_REF]. The conventional way to cope with this difficulty is to determine an optimal frequency band which maximizes the signal-tonoise ratio (SNR)i.e. the energy of the transients with respect to that of the background noisebefore enveloping (incidentally, this is also theoretically required to properly define the "envelope" of a signal [START_REF] Picinbono | On instantaneous amplitude and phase of signals[END_REF]). This issue has been well addressed and has led to the use of the SES in combination with different pre-processing tools apt to identify the optimal demodulation band, for instance rooted on the spectral correlation or, most frequently, on the spectral kurtosiscomputed with the fast kurtogram-as an effective measure of the 'impulsiveness' hidden in a signal [START_REF] Randall | Rolling element bearing diagnostics-a tutorial[END_REF][START_REF] Antoni | The spectral kurtosis: a useful tool for characterising non-stationary signals[END_REF][START_REF] Antoni | The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines[END_REF].The use of the kurtosis at this stage is further supported by the fact that it is related to the sum of the peaks in the SES divided by the zero-frequency SES [START_REF] Borghesani | The relationship between kurtosis-and envelope-based indexes for the diagnostic of rolling element bearings[END_REF]. A recent extension of the kurtogram, with similar goal but refined properties, is the infogram [START_REF] Antoni | The infogram: Entropic evidence of the signature of repetitive transients[END_REF][START_REF] Wang | An extension of the infograms to novel Bayesian inference for bearing fault feature identification[END_REF][START_REF] Li | Extracting repetitive transients for rotating machinery diagnosis using multiscale clustered grey infogram[END_REF]. Other approachespossibly used in conjunction with the latter ones-are based on first denoising the signals by various techniques before computing the SES [START_REF] Wang | Multiscale envelope manifold for enhanced fault diagnosis of rotating machines[END_REF][START_REF] Yu | Application of EMD method and Hilbert spectrum to the fault diagnosis of roller bearings[END_REF].

Although not strictly necessary for diagnosis, the final step consisting in extracting the fault signal is useful for assessing the severity of the fault and for designing health indicators for trending and for prognostics. From the signal processing point of view, this is yet a different task. One way to formulate the problem is as a blind deconvolution task, where the impulses due to the contact force on the defect are to be recovered from the measurements without knowledge of the impulse response [START_REF] Endo | Enhancement of autoregressive model based gear tooth fault detection technique by the use of minimum entropy deconvolution filter[END_REF][START_REF] Obuchowski | Blind equalization using combined skewness-kurtosis criterion for gearbox vibration enhancement[END_REF][START_REF] Mcdonald | Multipoint Optimal Minimum Entropy Deconvolution and Convolution Fix: Application to vibration fault detection[END_REF]. A different and probably less ambitious approach is to recover the transients (i.e. the responses to the impulses) as if they had been observed in the absence of background noise [START_REF] Boustany | Blind extraction of a cyclostationary signal using reduced-rank cyclic regression-A unifying approach[END_REF][START_REF] Cai | Sparsity-enabled signal decomposition using tunable Q-factor wavelet transform for fault feature extraction of gearbox[END_REF][START_REF] He | Sparsity-based algorithm for detecting faults in rotating machines[END_REF][START_REF] Zhang | Kurtosis based weighted sparse model with convex optimization technique for bearing fault diagnosis[END_REF][START_REF] He | Repetitive transients extraction algorithm for detecting bearing faults[END_REF][START_REF] Ding | Time-frequency manifold sparse reconstruction: A novel method for bearing fault feature extraction[END_REF].The former solutions actually require the use of sophisticated signal processing methods. A suboptimal solution justified by its simplicity is to extract an estimate of the fault signal by bandpass filtering with the optimal filter used in the computation of the SES. This brief tour of the state-of-the-art shows that bearing diagnosis often involve a combination of different signal processing methods which require a high level of expertise from the user. As a consequence, one of the current challenges is to develop a complete standalone diagnosis methodology that can be run without (or very limited) intervention of the user. This is referred to as automated diagnosis, or "semi-automated" to highlight the fact that even though the processing is fully automatic, the final decision is still taken by the end-user upon interpretation of the results. (Semi-) automated diagnosis of rolling element bearings has been addressed from two main points of view. The first one is by means of classification algorithms borrowed from the domains of data mining and machine learning. Although a large body of literature exists on the subject [START_REF] Srividya | Automated diagnosis of rolling element bearing defects using timedomain features and neural networks[END_REF][START_REF] Cococcioni | Automatic Diagnosis of Defects of Rolling Element Bearings Based on Computational Intelligence Techniques[END_REF][START_REF] Castejón | Automated diagnosis of rolling bearings using MRA and neural networks[END_REF][START_REF] Gryllias | Automated diagnostic approaches for deffective rolling element bearing using minimal training pattern classification methods[END_REF][START_REF] Lei | An intelligent fault diagnosis method using unsupervised feature learningtowards mechanical big data[END_REF][START_REF] Khan | Automated Bearing Fault Diagnosis Using 2D Analysis of Vibration Acceleration Signals under Variable Speed Conditions[END_REF][START_REF] Dibiano | Automated diagnostics for manufacturing machinery based on well-regularized deep neural networks, Integration[END_REF], it seems that the majority of the proposed solutions are confined to apply to a specific installationon which a classifier has been trained -and cannot be transposed to data recorded on other systems which have not been incorporated in the learning set. The second point of view proceeds from the signal processing side and ambitions a greater generality than the first one, yet it also appears much more challenging. There is indeed a very limited number of publications on the subject and they rarely address all the aspects of the diagnosis methodology [START_REF] Yang | Basis pursuit-based intelligent diagnosis of bearing faults[END_REF][START_REF] Rafiee | Application of mother wavelet functions for automatic gear and bearing fault diagnosis[END_REF][START_REF] Chen | A study on automatic machine condition monitoring and fault diagnosis for bearing and unbalanced rotor faults[END_REF][START_REF] Yang | Automatic Condition Monitoring of Industrial Rolling-Element Bearings Using Motor's Vibration and Current Analysis[END_REF][START_REF] Gerber | Time-Frequency Tracking of Spectral Structures Estimated by a Data-Driven Method[END_REF][START_REF] Firla | Automatic characteristic frequency association and allsideband demodulation for the detection of a bearing fault[END_REF][START_REF] Wang | Intelligent fault diagnosis of rolling element bearings using sparse wavelet energy based on overcomplete DWT and basis pursuit[END_REF].One methodology which has proved quite successful is described in Refs. [START_REF] Randall | Rolling element bearing diagnostics-a tutorial[END_REF][START_REF] Smith | Rolling element bearing diagnostics using the Case Western Reserve University data: A benchmark study[END_REF][START_REF] Sawalhi | Semi-automated bearing diagnostics-three case studies[END_REF] and has been further elaborated in Refs. [START_REF] Borghesani | Application of cepstrum pre-whitening for the diagnosis of bearing faults under variable speed conditions[END_REF][START_REF] Peeters | A comparison of cepstral editing methods as signal pre-processing techniques for vibration-based bearing fault detection[END_REF]. It consists in first pre-processing the signal, typically by whitening, in order to enhance the presence of low-energy transients. Next, an optimal band for demodulation is found from the fast kurtogram, from which the SES is computed. Finally, the fault signal is estimated by bandpass filtering in the optimal band returned by the kurtogram. In its trimmest version (where whitening is achieved by cepstral editing [START_REF] Borghesani | Application of cepstrum pre-whitening for the diagnosis of bearing faults under variable speed conditions[END_REF]), the whole process requires only one parameter to set which is the decomposition depth of the kurtogram. In many instances, the latter can actually be set by default to a predefined value. More sophisticated versions essentially improve the pre-processing step by using other whitening schemes [START_REF] Barbini | Phase editing as a signal pre-processing step for automated bearing fault detection[END_REF][START_REF] Peeters | International Congress on Sound & Vibration[END_REF] or blind deconvolution [START_REF] Sawalhi | The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis[END_REF]. Eventually, statistical tests can be designed on the so-obtained SES by following the lines of Refs. [START_REF] Antoni | Cyclic spectral analysis in practice[END_REF][START_REF] Girondin | Bearings fault detection in helicopters using frequency readjustment and cyclostationary analysis[END_REF][START_REF] Borghesani | Cyclostationary analysis with logarithmic variance stabilisation[END_REF][START_REF] Borghesani | CS2 analysis in presence of non-Gaussian background noise -Effect on traditional estimators and resilience of log-envelope indicators[END_REF] in order to automatically detect a fault with a given probability of false alarm.

The aim of the present work is to proceed with the objective of semi-automated diagnosis by using more advanced signal processing tools, while still simplifying the pre-processing step and keeping the number of tuning parameters to its minimum. A solution is proposed which addresses at once, in the same algorithm, the three goals of fault detection, fault identification and fault extraction. It is rooted on a shorttime-Fourier-transform (STFT) representation of the signal and therefore requires only one parameter setting corresponding the spectral resolution. As compared to the aforementioned methodology (whitening + fast kurtogram + SES + filtration), it has the advantage of not requiring any pre-processing (e.g. before computing the STFT), of being optimal for detectionin the Neyman-Pearson senseand of extracting a full-band version of the fault signal. It is also interesting as such since it achieves similar goals from a different way of processing the data and therefore offers a methodological diversity that is often desirable to improve the robustness of a diagnostic system. The use of the STFT is justified because it captures well the time-frequency structure of a series of transients produced by a bearing fault. Since the objective is not to arrive at a visual detection, there is no critical limitation to expect from the uncertainty principlee.g. as compared to other time-frequency representations, with better time-frequency localization, which have been advocated for diagnosis [START_REF] Feng | Recent advances in time-frequency analysis methods for machinery fault diagnosis: a review with application examples[END_REF] as long as the STFT returns the whole information contained in the analyzed signal, that is as long as it is invertible. In addition, the STFT is a linear transform that can be computed by means of very efficient algorithms and is therefore well suited to semi-automated diagnosis. The novelty of the present work is to endow the STFT coefficients of a bearing signal with a probabilistic model which switches between states: one where only background noise is observed and one where it is superposed with the occurrence of a transient. The label of the state is represented by a random latent variable which is encoded into a hidden Markov model (HMM). Since the STFT coefficients quickly tend in distribution to a complex-valued Gaussian in virtue of the Central Limit theorem [START_REF] Brillinger | Time Series: Data Analysis and Theory[END_REF], this is eventually formulated as a Gaussian mixture model. It thus happens that the value of the latent variable is related to a generalized likelihood ratio (LLR), from which a detection test can be designed that is optimal in the Neyman-Pearson sensei.e. which maximizes the probability of detection of a fault for a given probability of false alarm. Besides, the Fourier transform of the LLR provides a spectrum which is in all point similar to the SES, yet obtained automatically without need for careful pre-filtering. Finally, after the GMM has been identified, the corresponding faulty state can be extracted as a byproduct and the corresponding time signal recovered by inverting the STFT. It is highlighted that the fault signal is extracted in full-band, contrary to the output of the fast-kurtogram. Ultimately, it is shown that the proposed approach is general enough to deal with nonstationary operating conditions.

It is noteworthy that HMMs have been proposed in 1970s as statistical modelsfor time series and so far they have been applied in a wide range of fields including speech recognition, computer vision, pattern recognition and many other areas. Numerous works based on HMMs have been reported in fault diagnosis during the last decade [START_REF] Li | Hidden Markov model-based fault diagnostics method in speed-up and speed-down process for rotating machinery[END_REF][START_REF] Tobon-Mejia | A mixture of gaussians hidden markov model for failure diagnostic and prognostic[END_REF][START_REF] Miao | Condition monitoring and classification of rotating machinery using wavelets and hidden Markov models[END_REF][START_REF] Ocak | HMM-based fault detection and diagnosis scheme for rolling element bearings[END_REF][START_REF] Heyns | Combining synchronous averaging with a Gaussian mixture model novelty detection scheme for vibration-based condition monitoring of a gearbox[END_REF][START_REF] Ocak | A new bearing fault detection and diagnosis scheme based on hidden Markov modeling of vibration signals[END_REF][START_REF] Tobon-Mejia | A data-driven failure prognostics method based on mixture of Gaussians hidden Markov models[END_REF][START_REF] Yu | Health condition monitoring of machines based on hidden markovmodel and contribution analysis[END_REF][START_REF] Schmidt | A novelty detection diagnostic methodology for gearboxes operating under fluctuating operating conditions using probabilistic techniques[END_REF][START_REF] Zhou | Detection and diagnosis of bearing faults using shift-invariant dictionary learning and hidden Markov model[END_REF], however they are mainly concerned with the use of HMMs as classifiers [START_REF] Li | Hidden Markov model-based fault diagnostics method in speed-up and speed-down process for rotating machinery[END_REF][START_REF] Tobon-Mejia | A mixture of gaussians hidden markov model for failure diagnostic and prognostic[END_REF][START_REF] Miao | Condition monitoring and classification of rotating machinery using wavelets and hidden Markov models[END_REF][START_REF] Ocak | HMM-based fault detection and diagnosis scheme for rolling element bearings[END_REF][START_REF] Zhou | Detection and diagnosis of bearing faults using shift-invariant dictionary learning and hidden Markov model[END_REF][START_REF] Boutros | Detection and diagnosis of bearing and cutting tool faults using hidden Markov models[END_REF]and not for modelling the vibration signal itself as proposed in the present work.

The paper is organized as follows. Section 2 first introduces the probabilistic model based on the STFT decomposition and next explains the inference of the model parameters by means of the EM (expectedmaximization) algorithm. Section 3 provides the semi-automated diagnosis methodology that addresses the issues of optimal fault detection by designing a generalized likelihood ratio test (GLRT), of fault identification by means of the Fourier spectrum of the LLR and of fault extraction by means of a timevarying filter obtained from inverting the STFT. Section 4 addresses the important issue of parameter settings and of algorithm initialization which is hereafter verified by synthetic signals in section 5. Finally, section 6 validates the proposed methodology on several vibration signals and compare it with the reference methodology based on the (whitening + fast kurtogram + SES + filtration) sequence.

Probabilistic model

This section introduces the probabilistic model, its corresponding assumptions, and the inference of its parameters by the EM algorithm.

Signal model and STFT decomposition

Let denote the measured signal in the time domain, the part that contains the diagnostic information and the background noise. By definition, the "informative signal" and the "noise" are assumed mutually independent. The measured signal is thus expressed by the additive model .

(

) 1
Although the background noise actually comprises multiple components, it intervenes in model (1) as a global noise process, (e.g. with an equivalent covariance matrix). It is noteworthy that its exact probability distribution is not needed at this stage, yet a fair (and widely accepted) assumption is to model it as stationary. In contrast, since represents the fault signal in its early stage, it is well modelled by a series of impacts that repetitively excite resonances of the bearing and of its receiving structure, thus leading to successive damped impulse responses. Having a localized signature both in time and in frequency, such transients are well captured in a time-frequency decomposition, on the contrary to the stationary background noise which is spread all over the time-frequency plane. Although several time-frequency decompositions are possible, the proposed approach only requires an invertible one. The STFT meets this property while being associated with efficient algorithmic implementations. In addition, due to the Central Limit Theorem applied to the DFT (discrete Fourier Transform), the coefficients of the STFT quickly tend in distribution to a complex-valued Gaussian [START_REF] Brillinger | Time Series: Data Analysis and Theory[END_REF], an assumption that will substantially simplify the probabilistic model introduced hereafter. It is noteworthy that the aim is to decompose the signal without loss of information in the time-frequency domain in which it will be processed; therefore, considerations such as the uncertainty principle will not matter because they are relevant to time-frequency representation (i.e. visual analysis), which is not of concern here.

The STFT of signal over a time interval of length is defined as

(2)
where denotes a positive and smooth -long data-window which truncates a segment of the -long signal at time datum ( , ) with window shift R (1 <R< ) and where denotes the frequency (from 0 to ) with frequency resolution and bin index = 1, …, with .

It is noteworthy that is related to the "instantaneous complex envelope" of the signal described in both time and frequency. More precisely, its squared magnitude reflects the energy flow which is mapped by time index and frequency centered in a narrow frequency band [START_REF] Antoni | Cyclostationarity by examples[END_REF].

The next subsection introduces a two-state HMM to account for the different probability distributions of the STFT coefficients depending on whether a transient is present or not.

Hidden Markov model

Hereafter, the STFT coefficients are collected for all frequency bins at a given time instant i in a

Since the latent variable can take only two values, 0 or 1, with a given probability, say , it follows (by definition)the Bernoulli distribution, :

Here means the presence of noise only, i.e. "State 0: ", whereas indicates the presence of noise and the signal of interest, i.e. "State 1:

". This is recognized as a HMM with two hidden states.

Let us now introduce the likelihood function, ,

where ) denotes the circular-symmetric complex normal distribution with mean and covariance matrix which naturally arrives as a result of the Central Limit theorem applied to the STFT coefficients [START_REF] Brillinger | Time Series: Data Analysis and Theory[END_REF]. Without loss of generality, it is assumed that (as obtained after first centering the signal). Although the covariance matrix will be full in general, it has been observed in numerous cases that a diagonal structure provides an excellent approximation. Diagonal covariance matrices will therefore be assumed from now onwards, with the advantage of considerably limiting the number of unknowns in the model as well as providing a clear physical interpretation to the diagonal elements in terms of the instantaneous spectral density of signal y(t) (at time i).

Meanwhile, all the unknown parameters of the proposed model are denoted as . It is highlighted that the latent variables are hidden in the sense that they are not observed directly. Assuming independent segments in the STFT, the complete log-likelihood function is evaluated from Eq. (6) as .

Developing further, one has

where it has been assumed that all states are apriori equally probable.

The parameters are next estimated by maximizing the above likelihood function. In theory, this completely solves the problem since the estimated latent variable will then return the times of occurrence of the impacts on the faults and therefore the bearing characteristic frequency. Since it is difficult to find a closed-form solution, the EM algorithm [START_REF] Bilmes | A gentle tutorial of the EM algorithm and its application to parameter estimation for Gaussian mixture and hidden Markov models[END_REF] is used as an iterative method to find the maximum likelihood estimates. The EM algorithm makes use of the following quantities.

First, the posteriori probability distribution of the latent variable is formed as . [START_REF] Lai | Weak-signal detection based on the stochastic resonance of bistableDuffing oscillator and its application in incipient fault diagnosis[END_REF] According to Eq. ( 10), the expectation of the latent variable given the measurement is then computed as [START_REF] Guo | Multi-frequency weak signal detection based on multi-segment cascaded stochastic resonance for rolling bearings[END_REF] where denotes . After simple arrangement, this is expressed as [START_REF] Li | Gear fault detection using multi-scale morphological filters[END_REF] where denotes the natural logarithm of the likelihood ratio between the component of interest and the noise as given by .

(

) 13
The EM algorithm is summarized in Table 1.It is noted that it involves two parameters: for the maximum number of iterations and for the expected relative tolerance between and . These can be easily set by default. The estimation of the covariance of the signal of interest is finally obtained as where operator keeps only the positive eigenvalues of a matrix (here the positive elements of a diagonal matrix).

With these estimated parameters, an automatic fault detection scheme and a time-varying filter for filtering out the signal of interest are proposed in the next section.

b) K-component case

Let us now consider the general case with components of interest, . In this case, the probability that two (or more) components occur together will be assumed so small that such events will be disregarded(this may be seen as an extreme sparse representation where only one state is allowed at a time). Therefore, under the mutually exclusive assumption, the occurrence of the state is defined as [START_REF] Antoni | Cyclic spectral analysis in practice[END_REF] where indicates the presence of the signal . The pure noise casei.e. }is denoted as the state. Thus, there is a total of possible states in the model. Therefore, the marginal probability distribution treads [START_REF] Antoni | Fast computation of the spectral correlation[END_REF] where the latent variable follows the Bernoulli distribution, :

.

(Note that the latent variables are now dependent because of the mutually exclusive assumption.)

The posteriori probability distribution then reads

. ( 17
)
Assuming mutually exclusive states, the expectation of the latent variable is thus

. ( 18
)
Figure 1 illustrates the situation with K = 2 components, which involves three states.

Proposed methodology

This section now explains how the proposed HMM can be used for diagnosis. As explained in the introduction, the diagnosis tasks are divided into the detection, the identification and the extraction of a fault.

Detection

Being base on a Gaussian mixture model, the proposed probabilistic model easily lends itself to the formulation of a statistical detection test. This is formally stated by comparing two hypothesis, the null one, H 0 , corresponding to the situation where no transient is present in the signal and the alternative one, H 1 , where transients (and therefore a fault) are present. The test of hypothesis therefore literally reads [START_REF] Antoni | Cyclic spectral analysis of rolling-element bearing signals: facts and fictions[END_REF] or, in mathematical language,

An advantageous solution for testing H 1 against H 0 is by means of the generalized likelihood ratio test (GLRT). The principle is to take the ratio (or its logarithm) of the probability under H 1 to the probability under H 0where all unknown parameters are replaced by their maximum likelihood estimatesand to reject the null hypothesis if the ratio happens to be greater than a certain value. The GLRT has the advantage of being optimal in the Neyman-Pearson sense, that is it maximizes the probability of detection (here of accepting H 1 when there is actually a fault) for a given probability of false alarm (accepting H 1 when there is no fault). Another definite advantage of the GLRT in the present case is to come with a theoretical value for the threshold against which to compare the ratio.

Using Eq. ( 20), the logarithm of the GLRT is [START_REF] Abboud | Envelope analysis of rotating machine vibrations in variable speed conditions: A comprehensive treatment[END_REF] where , and are the maximum likelihood estimates returned by the EM algorithm (see Table 1) and [START_REF] Ming | Envelope calculation of the multi-component signal and its application to the deterministic component cancellation in bearing fault diagnosis[END_REF] Because the hypotheses H 0 and H 1 are nested, it can be shown from Wilk's theorem that, under H 0 , twice the logarithm of the GLRT follows a Chi2 distribution with number of degrees of freedom equal to the difference between the number of unknown parameters under H 1 (1 for , for and for ) and

H 0 ( for ), viz (23)
Therefore, the test of hypothesis is Reject H 0 if at the risk  (i.e. probability of alarm), where is the quantile of the Chi2 with probability 1-.

Identification

As explained in the introduction, the main tool to identify a bearing fault is the envelope spectrum due its ability to estimate the bearing characteristic frequencies. In the proposed HMM, the LLR in Eq. ( 13) may be interpreted as the likelihood that a transient is present at any time datum i and its variation in time reveals the fault frequency. Therefore, the Fourier transform of the LLR turns out a valid alternative to the SES. Other quantities also contain the information on the fault frequency such as the latent variables and the logarithm of the GLRT in Eq. ( 21), yet it has been verified in numerous examples that the spectrum of the LLR is more accurate than the former and just as good as the latter while involving a simpler expression. One shared advantage of the LLR and GLRT is that they can capture subtle modulations of the transient magnitudes (as would typically happen when the fault is subjected a nonuniform load distribution), which is less likely the case for the binary latent variable.

Fault extraction

As shown in this subsection, the fault signal can be reconstructed in full band based on the latent variables and the covariance matrices estimated in the HMM. First, let introduce the posterior probability distribution of the signal of interest as [START_REF] Antoni | The spectral kurtosis: a useful tool for characterising non-stationary signals[END_REF] The posterior probability density at a given frequency then reads [START_REF] Antoni | The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines[END_REF] where and stand for the noise and signal variance of the k th component frequency , respectively. After some manipulations, Eq. ( 25) is expressed as [START_REF] Borghesani | The relationship between kurtosis-and envelope-based indexes for the diagnostic of rolling element bearings[END_REF] with [START_REF] Antoni | The infogram: Entropic evidence of the signature of repetitive transients[END_REF] Therefore the expectation of the signal of interest is .

Finally, the time signal is obtained from Eq. ( 28) by using the inverse STFT.

Two remarks are noteworthy. First, it is seen that Eq. ( 28) corresponds to a time-varying filter from which superior performance is expected than from a conventional time-invariant filter. Second, the standard Wiener filter appears as a particular case under the assumption of stationarity, that is

, ( 29
)
where the latent variable for all time instants. In other words, Eq. ( 29) corresponds to the case where "State 1:

" occurs only.

The detection, identification and extraction tasks will be extensively illustrated in the next sections.

Parameter setting, algorithm initialization and validation

This section is divided in two parts. The first one addresses the setting of the parameters of the STFT decompositionwhich are the only ones not estimated in the proposed methodologyand in particular of the frequency resolution which is the only one that possibly requires an intervention by the user. The second part then addresses the initialization the EM algorithm by means of an effective data-driven approach.

Parameter setting

The parameters entering into the STFT are the window length and the window shift (see Eq. ( 2))

Window length

The value of directly controls the frequency resolution, ,

which characterizes the carrier frequency. It is required to cover at least the duration of a transient, which implies the condition

. ( 31
)
As the STFT is subjected to the uncertainty principle, , the highest switching frequency of the latent variable, , is bounded upward by [START_REF] Antoni | Fast computation of the kurtogram for the detection of transient faults[END_REF]. Therefore the available range of the latent variable is limited by

. ( 32
)
Therefore, should be taken short to allow a high switching rate in Eq. ( 32), but long enough to satisfy Eq. ( 31), i.e.

as illustrated in Fig. 2. Other reasons for taking small is to reduce the computation time required by the STFT and also to ensure sufficient segments for accurate parameter estimation.

It is noticed here that the rule for setting also works in the special case where the interval between adjacent transients is close to the transient duration . In other words, it is robust enough to balance the trade-off between a fine resolution and a high switching frequency of the latent variable. These facts will be further experimentally verified in subsection 4.3.

Window shift

There are two considerations for setting the window shift :

 first, for the STFT to be invertible, it is recommended to take at least overlap with a Hanning window,  second, if inevitability is not required, should be taken sufficiently small to keep enough diagnostic information while not increasing too much the computational cost and the dependence between adjacent segments; a typical choice is within and overlap with a Hanning window.

Therefore, the window shift can be easily set by default.

Initializing of the EM algorithm

The EM algorithm generally requires a good initialization for two reasons. The first ones is to avoid being trapped in possible local maxima of the likelihood and the second one is to achieve a fast convergence speed. A simple self-running solution is given hereafter to obtain good initial values for the covariance matrices of the two states in the HMM.

Initial estimates of the diagonal elements of the noise covariance matrix,

, are obtained by taking the median value of the natural logarithm of the squared magnitude of the STFT coefficients with respect to time instant , i.e.

. [START_REF] Obuchowski | Blind equalization using combined skewness-kurtosis criterion for gearbox vibration enhancement[END_REF] This estimator is based on the fact that State 0 is characterized by a higher probability than State 1 so that the median in the above equation is almost unaffected by the occurrences of the fault. Besides, the initial estimate of the covariance matrix in State 1 is based on the extreme values of .Specifically, let define the "deviation" [START_REF] Mcdonald | Multipoint Optimal Minimum Entropy Deconvolution and Convolution Fix: Application to vibration fault detection[END_REF] and a threshold equal to a fraction of the maximum of . Then, the sample set is constructed by collected all indices i where is found greater than . Hence, the diagonal elements of the covariance matrix in State 1 are estimated as [START_REF] Boustany | Blind extraction of a cyclostationary signal using reduced-rank cyclic regression-A unifying approach[END_REF] where is the cardinal of set . Finally, is obtained as .The corresponding probability is initialized to

. ( 36
)
It has been observed in numerous experiments that the proposed initializations are often quite close to the maximum likelihood estimates (global maximum)while allowing at the same time a fast convergence speed of the EM algorithm. This will be demonstrated in the next section.

Validation on synthetic signals

The section deals with a validation of the proposed methodology on synthetic signals.

5.1.Cases 1&2: demonstration of parameter selection

To demonstrate the performance of the proposed algorithm and its initialization, a synthetic signal is generated with a resonance frequency Hz, which is further modulated by a relatively high fault frequency Hz ( s, the sampling frequency is normalized as Hz). The synthetic signal is described as: (37) [START_REF] Zhang | Kurtosis based weighted sparse model with convex optimization technique for bearing fault diagnosis[END_REF] where and account for the uncertainties on the arrival time and on the magnitude of the transient, respectively. The white noise is set to a SNR of and the signal length is samples. A second-order system is defined by Eq. ( 38), whose numerator and denominator coefficients are and with , respectively. Figure3 shows the spectrogram (magnitude of the STFT) of the raw signal whose record in time is displayed in Fig. 4

(a).The frequency resolution is set to

Hz and a default value of 85% is used for overlap (R = 20). Following Eqs. ( 33)- [START_REF] Cai | Sparsity-enabled signal decomposition using tunable Q-factor wavelet transform for fault feature extraction of gearbox[END_REF], one can initialize the parameters , and as shown in Fig. 5 (a).It is seen that the proposed initialization is simple and effective, even though the estimated spectrum of the signal of interest still contains a significant contribution from noise especially below 0.08 Hz. After convergence of the EM algorithm, the estimation of the signal and noise spectra are close to the real values as can be seen in Fig. 5 (b). The very good estimation of the latent variable is verified in Fig. 4 (d). The detection capability of the HMM is now checked by means of the statistical test introduced in section 3.1. The same experience is repeated for several values of the signal-to-noise ratio (SNR) ranging from -20dB to 0dB and different values of the frequency resolutioncorresponding to . In each case, 500 hundreds realizations of the same signal (with different noise generations) are run. The signal length is 100,000 samples and the risk is set to 5%. The results are reported in Table 2. It is seen that the presence of the fault is detected for SNRs greater than or equal to -16dBfor all frequency resolution up to 2 -8 (SNRs even lower than this can be reached when the transient are lightly damped), which proves a certain robustness with respect to the latter parameter. It reduces to -8dB when the frequency resolution decrease to 2 -9 because it then becomes comparable to the fault frequency . For frequency resolution smaller than the fault could not be detected, which is consistent with the discussion of subsection 4.1.1. The results are also displayed in Fig. 6, which demonstrates that the proposed GLRT grows exponentially above the statistical threshold as the SNR increases. The statistical threshold is marked by a dotted horizontal line corresponding to a risk of 5%. The red bullet indicates the SNR where the GLRT goes above the threshold, i.e. where the fault is detected.

The same experiment is now carried on in order to check the effect of the initialization of the EM algorithm. The main parameter that governs initialization is the threshold  defined after Eq. ( 34). Different values are tested, i.e. , , , , which are represented in Fig. 7 (a) together with the deviation . ), which is probably the most difficult to obtain. It is seen that all initial estimates are very similar and also quite close the reference given by the squared magnitude of the frequency response . A closer look around the resonance in the band to actually shows that the initial estimation gradually improves when the threshold is increased.

The normalized root-mean-square error (RMSE) is displayed in Fig. 8 as a function of the iteration number k in the EM algorithm. It evidences that the convergence speed of EM algorithm also increases with the threshold. For information, thetotal computational time in the case is CPU (Central Processing Unit) with a PC with i7-3930K 3.20 GHz Processor. In order to verify the robustness with respect to the setting of the frequency resolution, a second synthetic signal is generated and analyzed with a deliberately coarse value of . The resonance frequency is Hz and the fault frequency Hz. The frequency resolution is chosen as ( Hz), which covers integer multiples of Hz but cannot resolve the resonance frequency . The corresponding limit on the detection of the fault frequency is Hz. All the other parameters are as in Case 1 (SNR = -6dB). Figure 9 shows the spectrogram of the raw signal. Since the modulation frequency is now close to the resonance one, this case encounters a trade-off between a fine spectral content and a large fault frequency range. After running the detection test, the GLRT is found equal to 3703 which is much greater than the statistical threshold at the risk of 5%,

. Hence the fault is detected without ambiguity. Figure 10 (a) displays the estimated LLR, which accurately localizes the fault occurrences: the function sharply takes very large values when it identifies a transient. Besides, the estimated latent variable locates exactly all the STFT segments that contains a fault occurrence, as shown in Fig. 10 (b). To further identify the fault type, the spectra of the LLR and of the latent variable are displayed in Fig. 11 together with the SES of the actual fault signal (note that all spectra are normalized to a unit maximum and that the limit is indicted by a vertical black dotted line). As explained in section 3, the spectrum of the LLR is preferred because it can capture modulations of the transient magnitudes which are not in the binary latent variable. It is indeed checked that the spectrum of the LLR in Fig. 11 (a) has a slightly larger extend than the spectrum of the latent variable in Fig. 11 (b). It is also noteworthy that the two spectra perfectly match the SES of the theoretical fault signal. Despite the coarse frequency resolution used in this case, this proves that the proposed method can still detect and identify the expected fault frequency with very good accuracy.

Case 3: Separation of a mixture of independent transients

This subsection demonstrates the potential of proposed HMM to deal with K = 2 simultaneous components in background noise. A synthetic signal with two different components is generated as shown in Fig. 12. The first component has a resonance frequency Hz with and its cyclic frequency is Hz ( samples), whereas the second one has a resonance frequency Hz with and cyclic frequency Hz ( samples)see Fig. 12 (a) and (b), respectively. All other parameters are set as in Case 1, with and . The SNR is 0 dB.

The spectrogram of the measurement computed with and (overlapping ratio ) is shown in Fig. 13. The estimated diagonals of the three covariance matrices are displayed in Fig. 14. It is seen that the spectra of all components are correctly identified. The corresponding latent variables and are displayed in Fig. 15. Compared with the reference signals, it is obvious that the times of occurrence of the two components have been correctly located. Finally, Fig. 16 displays the reconstructed repetitive transients and as well as their summation. Very good reconstruction is obtained, which demonstrates the performance of the proposed algorithm.

Validation on vibration signals

This section illustrates the application of the proposed methodology on actual vibration signals, where the three goals of detection, identification and fault extraction are again addressed separately. It is also demonstrated that the proposed methodology applies to time-varying operating conditions (data captured during a run-up), as often encountered in industrial applications. Comparison are made with the reference semi-automated diagnosis method described in the introduction where the signals are first whitened, then processed with the fast kurtogram in order to estimate an optimal frequency band for computing the SES and finally the fault is extracted by bandpass filtering.

Case 4: diagnosis of a ball fault

Three typical types of fault (i.e. inner race, outer race and ball fault) are investigated in a dataset from the Vibrations and Acoustics Laboratory of the University of New South Wales (Sydney) [START_REF] Antoni | Fast computation of the kurtogram for the detection of transient faults[END_REF]. The test-rig is a one-stage gearbox with primary and secondary shafts supported by ball bearings. Since it is often more difficult to identify a ball defect, particularly at incipient stage, this case is tested here. The spectrogram of the raw signal is displayed in Fig. 17. It is seen that there exists non-stationary components in the high frequency band above 10 kHz, whereas the low frequency range is dominated by high energy components related to the gearbox vibrations. The SES of the raw signal is displayed in Fig. 18. It is noted that there exists a relatively high value at the gearmesh frequency Hz surrounded by shaft speed sidebands ( Hz) which originates from the gearbox. Clearly, the information of the bearing fault is completely masked by high-energy components from the gearbox in the SES. At this point, it is therefore important to resort to methods that automatically select frequency bands in the signal where the SNR is maximum. As stated in Ref. [START_REF] Antoni | Fast computation of the kurtogram for the detection of transient faults[END_REF], the fast kurtogram has proved a powerful fourth-order spectral analysis tool for detecting and characterizing impulses in a signal. The fast kurtogram is applied here with decomposition levels in a -binary tree. As seen in Fig. 19, there exists several local maxima in the kurtogram. They are coherent with the spectrogram of Fig. 17 which evidences a clear non-stationary activity above 10 kHz. All the dyads with very high kurtosis values have been checked to have similar complex envelopes. Therefore one maximum is taken at 71.83 whose corresponding to the frequency band Hz. The SES in that band is displayed in Fig. 20. It clearly reveals the even BSF (ball spin frequency) surrounded by modulation sidebands at cage speed (FTF).

Next, the HMM model is estimated with the parameter settings listed in Table . The detection test returns a value of 171, 160 for the GLRT to be compared to a statistical threshold of 42 with a risk of 5%, which clearly concludes to the presence of a fault. Next, the LLR spectrum is displayed in Fig. 21. Comparing with Fig. 20, the LLR spectrum better enhances the odd harmonics of the BSF than the SES from the kurtogram, even if the diagnostics information is very similar in both cases. Before extracting the fault signal in the time domain, it is interesting to display the noise and signal spectra (diagonals of the covariance matrices and ). It is seen in Fig. 22 that the two spectra cross around 8 kHz, which is consistent with the two frequency bands identified in the spectrogram of Fig. 17. This reflects the fact that the high energy vibrations of the gearbox dominate the lower frequency range, whereas the repetitive transients dominate the higher frequency range. The capability of the proposed HMM method to reconstruct the full-band fault signal is now demonstrated and compared with the band-pass results obtained from the kurtogram.Figure23 (a) displays 1.36 s of the vibration signal and Fig. 24 (a) an enlarged view in the vicinity of a transient. The band-pass filtered signal in band Hz is displayed in Fig. 23 (b) and its enlarged view in Fig. 24 (b); it clearly evidences the presence of transients with maximum SNR. The reconstructed signal from the proposed HMM-based time-varying filter is displayed in Fig. 23 (c) and its enlarged view in Fig. 24 (c). As compared to the filtered signal based on the kurtogram, the reconstructed signal an exact location of the transients with their full-band spectral content. This may be used advantageously to better characterize the fault signature, infer the fault dimension and spectral content, and possibly update trend models for prognostics.

Cases5&6: diagnosis of bearing and gears

Signal recorded on another test rigare now considered. The test rig mainly consists of an electric asynchronous motor, a rotary encoder, 4 accelerometer sensors, a speed variator, a driving gear with 45 teeth, four bearings (3 healthy and 1 outer race fault) and two pinions (healthy and broken)see Fig. Two conditions are considered: 1) a damage bearing (outer race fault),and 2) a damage bearing (outer race fault) and a damaged gear together. It is highlighted here that the hole seen on the flank of the bearing in the upper right part of Fig. 25 is not the fault, but a through-bore used to seed a defect on the outer race. The actual size of the defect is about 2mm.

To demonstrate the effectiveness of the proposed method, two tests are considered hereafter: Case 5involvesa damage bearing (outer race fault) in cooperation with almost new gears and Case 6relates to the combination of the outer race fault and the broken pinion connected to the driving gear. The parameter settings for the HMM are listed in Table . The detection test returns a value of 392,330 for the GLRT to be compared to a statistical threshold of 42 with a risk of 5%, which clearly concludes to the presence of a fault. Figure 26 shows the spectrum of the LLR which indicates the harmonic structure of the suspected fault frequencies BPFO (Ball Pass Frequency on the Outer Race). Meanwhile the repetitive transients (characterized by the fault frequency),and the residual noise are separated from the raw signal as shown in Fig. 27.

b) Analysis of Case 6

Using the same test rig and the same parameter settings, a compound source of vibration that contains gear and bearing faults together is diagnosed in Case 6. The detection test returns a value of 117, 640 for the GLRT to be compared to a statistical threshold of 42 with a risk of 5%, which concludes to the presence of a fault. Figure28 shows the spectrum of the LLR which reveals the harmonics of BPFO and the harmonics of the secondary shaft ( ), thus demonstrating the simultaneous presence of the outer race defect and of the gear defect. The extracted transientssignal -and the residual noise -are displayed in Fig. 29 and their respective SES's in Fig. 30. Since the bearing fault has a transient nature it is recovered in whereas the gear fault has a more stationary nature (with slight modulations) that is recovered in .

6.3.Case 7: diagnosis of bearing in non-stationary operating conditions

This subsection illustrates the potential of the proposed methodology in variable operating conditions. A runup of from to has been manually produced with the test rig of Fig. 25. The instantaneous speed is displayed in Fig. 31 (a) and the corresponding acceleration signal in Fig. 31 (b). The latter undergoes speed-dependent variations in magnitude and phase. Since it places no constraint on the distribution of the time instants of the impacts, the HMM can in theory cope with machine signals recorded in time-varying operating conditions. However, in order to account for the speed-dependent variation of the probability distribution of the states, the whole signal is divided into consecutive speed segments of 12 rotations as suggested in Ref. [START_REF] Abboud | Speed-spectral whitening for enhancing envelope analysis in speed varying conditions[END_REF]-see Fig. 31 (c). Then, the parameters of the HMM are estimated separately on each segment. The probability of State 1 in the five operating conditions takes the following estimates: 0.28, 0.30, 0.30, 0.33 and 0.37.

The detection test returns values of consecutive speed segments, i.e. 749,520, 195,640, 161,570, 155,220 and 339,850 for the GLRT to be compared to a statistical threshold of 42 with a risk of 5%, which concludes to the presence of a fault.

Figs. 32-33 display the raw signal on a few selected segments together with the corresponding reconstructed transients. It is seen that the transients are all well identified, although their behavior slightly changes with the rotation speed: the 4 th segment corresponding to a speed around 20Hz exhibits the cleanest signature of the fault. The reconstructed components are displayed in Fig. 39 and 40 in intervals . While the bearing fault is identified (f 3 = BPFO) in the second component , the first component

Conclusion

The HMM provides a versatile tool for characterizing the non-stationarity of a signal from its secondorder statistics only, regardless of operating speed.

This paper has introduced a new stochastic model for representing the vibration signature of bearing faults in their incipient stage, when they have the typical form of a series of transients overlaid with stationary background noise. It is based on a HMM with intermittent switches between states characterized by Gaussian distributions with different covariance matrices. The main goal of the model is to provide a semi-automated diagnosis method that can handle at the same time the tasks of fault detection, fault identification and fault extraction. As compared to other semi-automated diagnostic solutions, the proposed one provides several advantages: first, the detection test is optimal in the Neyman-Pearson sense, second the spectrum of the LLRa quantity equivalent to the SES-does not require preprocessing, and third the extraction of the fault signal is full-band and based on a time-varying filter superior to alternative stationary filters. The performance of the proposed method has been demonstrated on several vibration signals from test rigs (gear and bearing). Superior or equivalent results to the conventional methodology based on the fast kurtogram and envelope analysis and have been observed in all cases. Results are clearly superior in terms of reconstruction of the fault signal when compared to conventional band-pass filtering. The proposed model also deals with the case where there exists multiple components and with signals acquired under non-stationary operating conditions of the machine. The proposed methodology requires the setting of only one parameter, the frequency resolution; it has been verified to be quite robust in this respect provided that the frequency resolution is set greater than the potential fault frequencies. The performance is however based on the validity of the stochastic model and there might be cases which cannot be represented by the latter; in particular, a marked limitation, is that the proposed HMM applies only to incipient/localized faults with a transient nature. Nevertheless, the methodology introduced in the paper may open the way to the proposal of other types of stochastic models that can capture fault configurations not covered in the present work.

  a noisy component . All components in the model are allowed to have different probability distributions and are controlled by a vector of latent variables, , each of which acting as a switch taking only values 0 or 1. One component case Let us first consider the simplest model with only one component of interest, that is

Figure 1 .

 1 Figure 1. Graphical model in the case of K = 2 components.

Figure 2 :

 2 Figure 2: Illustration of how to select the window length and shift with respect to transient durations T I and cycle t.

Figure 3 .

 3 Figure 3. Spectrogram of the signal simulated in Case 1 with resonance frequency Hz, and fault frequency Hz ( s, , and SNR = -6dB).

Figure 4 .

 4 Figure 4. (a) Synthetic signal of Case 1with white noise (SNR = -6 dB). (b) Synthetic repetitive transients. (c) Deviation (black line) with threshold (red dotted line) set to with .(d) Estimated latent variable .

Figure 5 .

 5 Figure 5. (a) Initialized spectra (diagonals of covariance matrices and ) (red line and blue asterisks) and (b) estimated ones from the EM algorithm (diagonals of covariance matrices and ) (red lines) together with the squared magnitude of the frequency response (black dashed line and blue asterisks).

Figure 6 .

 6 Figure 6. Evolution of the logarithm of the GLRT as a function of the SNR and the frequency resolution f = 1/N w .The statistical threshold is marked by a dotted horizontal line corresponding to a risk of 5%. The red bullet indicates the SNR where the GLRT goes above the threshold, i.e. where the fault is detected.

Figure 7 (

 7 b) displays the corresponding estimates of the transient spectrum in State 1 (i.e. the diagonal of covariance matrix

Figure 7 .

 7 Figure 7. (a) Deviation (black line) with threshold (colored dotted lines) set to , , , and (b) corresponding initialized spectrum (diagonal of covariance matrix) in colored lines together with the squared magnitude of the frequency response (black dashed line and blue asteriks).

Figure 8 .

 8 Figure 8. Normalized root-mean-square error (RMSE) on as a function of the iteration number of the EM algorithm ( and on ).

Figure 9 .

 9 Figure 9. Spectrogram of the signal simulated in Case 2 with resonance frequency Hz, and fault frequency Hz (SNR = -6dB).

Figure 10 (

 10 Figure 10(a) LLR(i)with the statistical threshold at the risk of 5%, and (b) latent variable in the time domain.

Figure 11 :

 11 Figure 11: Spectra of the LLR(blue line) (a) and of the latent variable (b) superimposed with the SES of actual fault signal (red dotted line). The limit Hz is indicted by a vertical black dotted line(Normalization to unit maximum value).

Figure 12

 12 Figure 12(a) Component with period samples ( and ).(b) Component with period samples ( and ).(c) Noisy measurement (SNR = dB).

Figure 13 .

 13 Figure 13. Spectrogram of the signal simulated in Case 3 with two components ( samples and samples respectively).

Figure 14 .

 14 Figure 14. Estimated spectra of component (red solid line), component (green solid line) and noise (blue dashed line). The squared magnitudes of the two frequency responses are indicated by black dashed lines.

Figure 15 .

 15 Figure 15. (a) Estimated latent variable (red solid line) together with the actual first component (black solid line). (b) Estimated latent variable (green solid line) together with the actual second component (black solid line).

Figure 16 .

 16 Figure 16. Reconstructed repetitive transients corresponding to (a) the first component with period samples and (b) the second component with period samples. (c) Summation of two components and .

Figure 17 .

 17 Figure 17. Spectrogram (logarithmic scale) of signal of Case 4 (frequency resolution Hz).

Figure 18 :

 18 Figure 18: Squared envelope spectrum of the raw signal (normalized to unit maximum value).

Figure 19 .

 19 Figure 19. Kurtogram of signal of Case 4 computed over levels with a -binary tree and a 8 coefficient prototype filter. Several local maxima are presented. One maximum is taken to the frequency band Hz.

Figure 20 .

 20 Figure 20. Squared envelope spectrum in frequency band Hz returned by the kurtogram (normalization to unit maximum value).

Figure 21 :

 21 Figure 21: Spectrum of the LLR with markers at the theoretical fault frequency and harmonics; the limit Hz is indicted by a vertical black dotted line(normalization to unit maximum value).

Figure 22 .

 22 Figure 22. Diagonals of covariance matrices (frequency resolutionHz); the spectrum of the fault signal is indicated by a red solid line and that of the noise by a blue dashed line.

Figure 23 (

 23 Figure 23(a) Vibration signal of Case 4 and (b) band-pass filtered signal in the frequency band Hz. (c) Full-band reconstructed fault signal from the proposed HMM-based time-varying filter.

Figure 24 (

 24 Figure 24(a) Enlarged view of vibration signal of Case 4 (with indication of the period corresponding to the peak in Fig. 18) and (b) enlarged view of band-pass filtered signal in the frequency band Hz. (c) Enlarged view of the full-band reconstructed fault signal from the proposed HMM-based time-varying filter.

25

 25

  .

Figure 25 .

 25 Figure 25. Test rig setup.

  Analysis of Case 5

Figure 26 :

 26 Figure 26: Spectrum of the LLR with markers at the suspected fault frequencies BPFOand its harmonics (normalization to unit maximum value).

Figure 27 (

 27 Figure 27(a) Vibration signal of Case 6 divided into (b)the full-band reconstructed fault signal from the proposed HMM-based time-varying filter and (c) the noise (residual) signal ().

Figure 28 :

 28 Figure 28: Spectrum of the LLR with markers at the suspected fault frequencies BPFO, its harmonics and (normalization to unit maximum value).

Figure 29 (

 29 Figure 29 (a) Vibration signal of Case 6divided into (b) the full-band reconstructed fault signal from the proposed HMM-based time-varying filter and (c) the noise (residual) signal ().

Figure 30 .

 30 Figure 30. Squared envelope spectrum of (a) the vibration signal of Case 6, (b) the full-band recovered time signal and (c) the noise (residual) signal (normalization to unit maximum value).

Figure 31 .

 31 Figure 31.(a) Estimated instantaneous speed of signal in Case 7 and (b) its corresponding acceleration signal which undergoes speed-dependent magnitude modulation.(c) Division of the estimated instantaneous speed in 5 operating conditions.

Figure 32 .

 32 Figure 32. (a) Raw signal in operating condition No. 2 and (b) the corresponding reconstructed transients from the HMM-based time-varying filter; (c) raw signal in operating condition No. 3 and (d) the corresponding reconstructed transients from the HMM-based time-varying filter.

Figure 33 . 3 Figure 35 .

 33335 Figure 33. (a) Raw signal in operating condition No. 4 and (b) the corresponding reconstructed transients from the HMM-based time-varying filter; (c) raw signal in operating condition No. 5 and (d) the corresponding reconstructed transients from the HMM-based time-varying filter.

Figure 36 .

 36 Figure 36. Spectrogram (logarithmic scale) of signal of Case 8 with evidence of two states: (a) global distribution with spectral content in band and (b) local distribution with spectral content in band and .

Figure 37 .

 37 Figure 37. Diagonals of covariance matrices of the three components, (red solid line), (green solid line) and noise (blue dashed line).

Figure 38 .

 38 Figure 38. Spectra of the LLR of the two components: (a) and (b) , respectively (normalization to unit maximum value).

Table 1 :

 1 Explicit steps of the EM algorithm to infer the parameters in the HMM.

Table 2 :

 2 parameters used in the detection test.

		2 -5 = 1/32	2 -6 = 1/64	2 -7 = 1/128	2 -8 = 1/256	2 -9 = 1/512
		15	31	63	127	255
	N (STFT length)	19994	9994	4994	2558	1293
	Threshold	13	23	42	78	147
	SNR from which H 0 is rejected	-16dB	-16dB	-16dB	-16dB	-8dB

Table 3 :

 3 Parameter settings in Case 4.

	Sampling frequency (Hz)	48000
	Duration (s)		1.365
			20
	Rotation frequency -	(Hz)	10
	Ball spin frequency -BSF (Hz)	26.11
	Fundamental train frequency -FTF (Hz)	4.08

Table 4 :

 4 Parameter settings in Case 5 and Case 6.

Acknowledgments

This work was supported financially by the Fundamental Research Funds for the Central Universities (Grant no. 2018RC009), the National Key Research and Development Program of China (Grant no. 2016YFB1200402), a scholarship from the China Scholarship Council (Grant no. 201304490040) and was performed within the framework of the Labex CeLyA of University of Lyon, operated by the French National Research Agency.

The identification of the fault characteristic frequency requires specific processing under non-stationary operating conditions. Since the defect impacts occur periodically with respect to the angular position, its frequency is to be computed in the order domain [START_REF] Abboud | The spectral analysis of cyclo-non-stationary signals[END_REF]. Therefore, the phase-corrected STFT in Eq. ( 2) have been resampled (using cubic splines interpolation) from the time to the angular domain, while maintaining a constant spectral bandwidth. The order spectrum of the LLR has then been computed on the resampled data. As seen in Fig. 34, it clearly reveals the presence of the Ball Pass Order on the Outer race (BPOO).

Case 8: diagnosis in the presence of multiple components

In order to demonstrate the performance of the HMM in the case of multiple-components, the dataset supported by the Department of Mechanical Engineering of Curtin University (Bentley) and made available online(http://data-acoustics.com/measurements/bearing-faults/bearing-1/) has been used. It corresponds to radial vibration measurements taken on the bearing housing of the SpectraQuest Machinery Fault Simulator test rig with a known outer race bearing fault. This case is interesting since there exists two different probabilities of states, as shown in Fig. 35. Their corresponding spectrograms are displayed in Fig. 36. Sampling frequency (Hz) 51200 Duration (s) 10 corresponds to transients that occur in the second time interval only; since their frequency happens to be about twice the cage speed, it might indicate that the bearing is also misaligned.