Multi-Band Frequency Window for Time-Frequency Fault Diagnosis of Induction Machines

Jordi Burriel-Valencia, Ruben Puche-Panadero, Javier Martinez-Roman, Angel Sapena-Bano*, Martin Riera-Guasp, Manuel Pineda-Sanchez The authors are with the Institute for Energy Engineering, Universitat Politècnica de València, Camino de Vera s/n, 46022, Valencia, Spain; jorburva@die.upv.es, rupucpa@die.upv.es, jmroman@die.upv.es, asapena@die.upv.es, mriera@die.upv.es, mpineda@die.upv.es * Correspondence: asapena@die.upv.es; Tel.: +34-96-3877597 † All the authors contributed equally to this work.


Introduction
Induction machines (IMs) are a key component of many industrial processes, either as motors or as generators, such as double fed induction generators (DFIGs) used for wind energy generation.
Their reliability ensures the continuity of the production processes, but they are subjected to eventual failures (broken bars, bearing faults, eccentricity, turn to turn or phase to ground short-circuits, etc.), which may cause unexpected breakdowns and high economic losses.A way to reduce these risks is the continuous monitoring of the machine's condition, in order to detect the presence of a fault at an early stage, when corrective measures can be implemented or maintenance works scheduled.
Bearing (inner race) Bearing (balls) A diagnostic technique that has gained a widespread interest in recent years is the motor current signature analysis (MCSA) [23][24][25][26], which is based on the detection of the characteristic fault signatures that each type of fault impress in the current spectrum.It is a non-invasive method, which can be applied on line without disturbing the normal operation of the machine; it requires, in its more basic implementation, just a current sensor for acquiring the current signal and a fast Fourier transform (FFT) for generating its spectrum; and it can detect a wide variety of machine faults through their spectral signatures, because each type of fault generate a specific set of fault frequencies, as shown in Table 1 [27], where f 1 stands for the fundamental component frequency, p is the number of pole pairs, s is the machine slip, f r is the rotational frequency of the IM rotor, N b is the number of balls of the bearings, D b if the bearing diameter, D c is the pitch or cage diameter, and β is the contact angle.
Nevertheless, the use of the current spectrum as signal processing tool in MCSA limits its field of application to machines working in stationary conditions, which is not the case of industrial processes with varying load or speed conditions, or also of wind generators operating under variable wind regimes.In these cases, the fault frequencies shown in Table 1 become random time functions, and the fault harmonics do not produce isolated spectral lines in the current spectrum, which blurs their characteristic signatures.
To extend MCSA to the fault diagnosis of IM working in transient regime, advanced time-frequency (TF) transforms of the current are needed, so that the transient fault signatures can be identified in a joint TF domain.These transforms can be lineal, such as the short-time Fourier transform (STFT) [28][29][30], the short-frequency Fourier transform (SFFT) [31], and the wavelet transform (WT) [32], or quadratic, such as the Wigner-Ville distribution (WVD) [33], or the ambiguity function [34].Quadratic TF transforms can achieve optimal resolution for mono-component chirp signals, but in case of multi-component ones, they produce cross-terms artifacts that pollute the TF representation of the current, making it difficult the correct identification of the fault harmonics.On the contrary, linear TF transforms are free from cross-terms artifacts.The STFT representation, the spectrogram, and the WVD representation, are built by multiplying the current signal with an analysing window at each time instant, and performing the FT of the resultant signal.In the case of the STFT, this window has a constant shape, which makes it difficult to obtain a good resolution both in time and in frequency: a short window gives a good time resolution, but a poor frequency resolution; on the contrary, a long window gives a good frequency resolution, but a poor time resolution.The WT solves this issue by performing a multi-resolution analysis, using different windows at different frequency bands: long windows for the lower frequency bands, and short windows for the higher frequency bands.But this approach makes the application of the WT much more complex than the FT, and distorts the TF signature of the fault harmonics.Diverse solutions to this problem have been proposed recently, in order to built a current's spectrogram with enough resolution for being used as an IM fault diagnosis tool.One of them is to adapt the shape of the analyzing window in the TF domain to the expected shape of the fault harmonics, either using a single window, as in [35], or using different window shapes for different sections of the current signal, as in [36].This approach requires a deep a priori knowledge of the time evolution of the fault harmonics in the current signal, which requires highly skilled maintenance personnel for implementing it, and difficult its application in automated diagnostic systems.Other approach apply to the spectrogram a post-processing based on reassignment [37] or synchrosqueezing [9,38] techniques, so improving its sharpness, but these techniques add a considerable computational burden to the process of building the current spectrogram, departing from the simplicity of the STFT.
Another alternative for obtaining an improved spectrogram is the Matching Pursuit approach [39].This method is based on calculating a series of spectrograms, using a set of different windows designated as "dictionary" [40], which has to be previously built.Then, combining spectrograms corresponding to each window of the dictionary through a pre-defined algorithm, obtains the final spectrogram, which is considered the optimum one.This method has several drawbacks, such as the need of building an extensive dictionary, with a great amount of different windows for obtaining a good resolution spectrogram.This implies to calculate a huge quantity of spectrograms and thus consuming a vast quantity of time and computational resources In this paper, a novel approach (up to the best of the authors' knowledge) is proposed to obtain a high-resolution current spectrogram, useful for fault diagnostic purposes, with the simplicity of a single STFT.It is based on 1. Performing the STFT with a wide range of windows with different lengths, and selecting, for each point in the TF domain, the best result obtained at that point among the complete set of windows.
2. Instead of running a separate STFT for each of the windows used in the analysis, a single, multi-band window is built by stacking all the desired analysing windows in consecutive frequency bands.This approach obtains in parallel the spectrograms corresponding to several hundreds of different analysis windows with the computing cost of a single one, which makes it suitable for fast, online diagnostic systems in transient regime.
The structure of the paper is as follows.In Section 2 the generation of the spectrogram with Gaussian windows shifted in the frequency domain is analyzed, and in Section 3 it is used for the theoretical and practical explanation of the proposed method.In Section 4, it is validated with the analysis of the start-up current of a high rated power, medium voltage squirrel cage induction motor, with broken bars.Section 5 presents the conclusions of this work.

Time-Frequency Analysis of the Machine's Current via STFT with a Multi-Band Window
To highlight the spectral content of a time-varying current signal, which may contain the characteristic fault harmonics given in Table 1, it is necessary to generate a representation of the current signal in the join TF domain.Among the diverse transforms available to this end (STFT, SFFT, WT, WVD, ambiguity function, etc.), the STFT has been selected in this work, because it is a linear transform, without cross-terms artifacts, and is computationally very effective, which makes it suitable for being implemented in low-cost, low-power embedded devices for on-line CBM systems.
In this Section, the traditional STFT analysis of the stator current using a Gaussian window will be first reviewed, and after the proposed method using a multi-band frequency window will be presented and compared with the traditional one.

Spectrogram of Machine's Current
The STFT of a current signal i(t) is a linear TF transform that is able to generate a joint TF representation of the current, the spectrogram, through the following steps: that emphasizes the content of the current signal at time τ, and attenuates it at other times 2. The Fourier transform if applied to the time-windowed signal i τ (t), which gives the frequency content of the current signal i(t) around time τ where ω = 2π f , and f stands for the frequency, in Hz.
3. The energy density spectrum at time τ is obtained as For each instant τ the STFT generates a different energy spectral density I SP (τ, ω), and the total set of these spectra constitutes the current spectrogram.A critical issue for obtaining a high resolution spectrogram of the current signal is the selection of the window h(t) in (4).To obtain a high resolution of the energy content of the current signal in the joint TF domain, it is necessary to use a window with a high concentration of energy in the TF plane, but such energy concentration is limited by the Heisenberg's uncertainty principle : a short time window gives a good time resolution, but a poor frequency resolution, and, on the contrary, a long time window gives a good frequency resolution, but a poor time resolution.The window that can achieve the highest energy concentration in the joint TF domain is the Gaussian window [35], given in the time domain by and in the frequency domain by The standard deviation of the Gaussian window in the time domain ( 5) is σ 2 t = 1/(2α), and in the frequency domain ( 6) is σ 2 ω = α/2.Therefore, for the Gaussian window, the product of its duration σ t and its bandwidth σ ω gives [41] The Heisenberg's uncertainty principle states that one cannot construct any signal whose duration σ t and bandwidth σ ω are, both, arbitrarily small, because In this way, the Gaussian window has a duration-bandwidth product (7) that reaches the minimum value (i.e, the highest concentration in the joint TF plane) that can be achieved under the uncertainty principle (8).
The parameter α in ( 5) and in ( 6) is the only one that defines the shape of the Gaussian window.A low value of α gives a long window with a narrow bandwidth, while a high value of α gives a short window, with a wide bandwidth.This parameter must be tuned to the current signal to be analyzed with the STFT.As detailed in [35], the optimal Gaussian window to build the spectrogram of a given current signal is the one that has the maximum overlap with the current in the TF domain.That is, the optimal parameter α is the one whose height/width ratio σ ω /σ t = α best approximates the slope of the current signal in the TF domain.Unfortunately, the slope of the fault harmonics in an IM in transient regime is not a constant value, because in this regime the slip and the rotational frequency in Table 1 are time-varying quantities, as well as the fault frequencies that depend on them.Besides, different components of the current signal may have different slopes (such as the fundamental component and the fault harmonics).These facts preclude the use of a single, optimal Gaussian windows for building a high resolution diagnostic spectrogram of the IM current.

Frequency Shifting of the Gaussian Analysing Window
If the Gaussian window ( 5) is used as the analysing window to build the energy density spectrum, then (4) becomes The Gaussian window ( 5) is a real valued function.If it shifted in the frequency domain by a frequency f k , corresponding to an angular frequency ω k = 2π f k , then (5) becomes a complex-valued function Replacing (10) in (9) gives that is, and, making the change of variable ω = ω + ω k , ( 12) can be expressed as Therefore, as the shift ω k of the Gaussian analysis window changes, the frequencies of the entire corresponding spectrogram (13) suffer the same ω k shift, at each time τ.

Proposed Multi-Band Analysing Window
A possible approach for obtaining a high-resolution spectrogram of a current signals with components of different slopes in the TF domain could be a variant of the Matching Pursuit approach, using as dictionary a set of Gaussian functions with different values of α.In this way, a batch of spectrograms, each one with a different value α for the Gaussian window (5) would be generated.
In a second stage, for each point in the TF domain, the best value obtained among the whole set of spectrograms would be selected, giving the best approximation to the ideal TF representation of the current signal.
The drawback of this technique is the high amount of resources that it requires.For each possible value of α, in a given range, a full spectrogram must be built, which is a time-consuming operation, and, also, it must be stored, which implies high memory resources.Afterwards, a processing algorithm must be applied to each point of all the spectrograms to combine them.These requirements make this technique unsuitable for being deployed with on-line, low power embedded devices.
On the contrary, the novel technique presented in this paper can achieve the same results at What is proposed in this paper is to extend this feature, using, instead of two frequency bands, a partition of the whole spectrogram in N g adjacent bands, each one with a frequency width equal to f b ; and to use a set of N g Gaussian windows, with different values of α (α k , k = 0 . . .N g − 1), and with increasing shifting frequencies ( , for filling each of these frequency bands.This can be achieved in a single run of the FFT if, using the superposition principle, all these N g , frequency shifted, Gaussian windows are summed in the time domain, giving a single time window. The use of this new, multi-band window, in (4) gives a single spectrogram with N g adjacent bands, each one corresponding to the analysis of the current signal with a different parameter α k for the Gaussian window.
For a given sampling frequency f s , and a frequency band of interest f b , the total number N g of adjacent frequency bands that can be used is equal to The proposed multi-band window is built, using (10), as As an example, a multi-band window (15) has been built using N g = 10 different Gaussian windows with values of α in (5) spanning linearly the range [α min = 1α max = 700], with a step equal to ∆α = (α max − α min )/(N g − 1), as The multi-band window ( 16) is displayed in Fig. 1, which shows the real (Fig. 1.top) and the imaginary part (Fig. 1.middle) of the window, as well as its spectrum (Fig. 1.bottom).Fig. 2 shows the spectrogram of the window in the joint TF domain, with all the individual Gaussian atoms stacked in adjacent frequency bands.

Steps for Applying the Proposed Multi-Band Window
In this sub-section the process for applying the concept of multi-band window to a given current signal i(t) is explained, using a synthetic current signal i(t) with a sinusoidal component of 50 Hz, and a linear chirp with a frequency slope of 1 Hz/s, starting at -50 Hz.It has been built with a duration of 50 s, using a frequency sampling rate of 200 Hz (Fig. 3), i(t) = cos(2π50t) + 0.05 cos(2π(−50t + t 2 )).The steps for analysing the current signal ( 17) using the proposed multi-band windows are the following ones: 1.The analysis window ( 15) is built: which gives the maximum number of elementary Gaussian windows from ( 14) (Ng = 200/100 = 2).The current signal is low-pass filtered with a cut-off frequency equal to f b .In this work, a spectral filter which zeroes all the spectrum bins with a frequency greater than f b has been used, as in [42].
• second, the parameters α k of each of these windows in (15) must be selected.For the simple signal (17), only two values of α are used, α 0 = 1 (a long window), and α 1 = 12.6, tailored to the chirp component according to [35].
The resulting multi-band window expression, applying (15), is This window is plotted in Fig. 4, in the time and frequency domains.

Experimental Validation
For the experimental verification of the proposed approach, an IM whose characteristics are given in Appendix A has been prepared with a forced rotor fault, by drilling a hole in one of the rotor bars.
The stator current during a startup transient has been acquired with a frequency rate f s = 5000 Hz, during 12 seconds, using a current probe whose characteristics are given in Appendix B, and it is shown in Fig. 8. From Table 1, the characteristic frequencies of the fault harmonics in an IM with a rotor asymmetry, such as a broken bar, are, for the main fault harmonics (k/p = 1 in Table 1), In particular, the fault harmonic with a frequency given by which is known as the lower side-band harmonic (LSH), is commonly tracked for the diagnosis of rotor asymmetries.During a start-up transient, the trajectory of the LSH in the TF plane given by (20) generates a typical V-shaped fault signature [43], with a frequency that initially (s = 1) is equal to the fundamental frequency f 1 , decreases to 0 (s = 0.5), and then increases again until its steady-state regime value of f 1 (1 − 2s) ≈ f 1 , according to (20).The ability to detect this fault harmonic using the proposed method is to be assessed in this section.
The three steps presented in Section 3.1 will be followed in this experimental case to obtain the spectrogram of the current signal presented in Fig. 8.
1.The analysis window ( 15) is built: • first, the bandwidth of diagnostic interest is established.In this case, from (20), the frequency band of interest is the [0 -50 Hz] band.Due to the presence of higher order harmonics in the current spectrum, apart form the LSH, a wider band [0 -125 Hz] has been selected, in order to better assess the strength of the LSH compared with them.In this way, With a sampling frequency of 5000 Hz, the maximum number of elementary Gaussian windows from ( 14) is Ng = 5000/125 = 40 windows.
• second, the parameter α k of each of these windows in ( 15) is selected.In this case, a linear range of the parameter α is used, covering the range [α min = 1α max = 700].
The resulting multi-band window, applying (15), is with The multi-band window (21) is displayed in Fig. 9, which shows the real (Fig. 9.top) and the imaginary part (Fig. 9.middle) of the window, as well as its spectrum (Fig. 9.bottom).

Conclusions
In this work, a novel technique for performing the fault diagnosis of IMs working in transient regime has been presented and validated experimentally.This technique consists in building a multi-band analysing window, composed by several, different Gaussian windows, stacked in the frequency domain.When this window is multiplied by the current signal in the time domain, it generates a spectrum which contains all the spectra generated by each of the Gaussian windows that form the multi-band window.In this way, a spectrogram containing even hundreds of different Gaussian windows can be obtained at the cost of a single run of the STFT algorithm.The selection of the parameters of the individual windows that compose the multi-band window can be setup by the user using different criteria, without affecting the performance of the proposed approach.In this work, a blind approach has been used, by selecting a range of individual windows that cover a wide range of, a priori, unkown fault harmonic components.This approach can be particularly useful for automated diagnostic system, which can operate without expert assistance, and for the detection of different types of faults.In this case, a single multi-band window can avoid the design of multiple analysis windows specifically designed for each single type fault.
roughly the cost of a single STFT, in terms of speed and storage requirements, even with the use of several hundreds of Gaussian windows with different values of α.It is based on a particular feature of the IM faults presented in Table1: in most industrial IMs, the fault harmonics with the highest amplitudes are those with an index k = 1, and they are located in a narrow, low frequency band with a bandwidth f b (normally of one or two hundreds of Hz).Nevertheless, the current signal is acquired normally using high frequencies rates, from 5 or 10 kHz up to 100 kHz and more.This implies that, when performing the FT of the windowed current signal, at each stage of the STFT process, only the values within the narrow band [0 − f b ] of diagnostic interest are kept, and the rest of the spectrogram values are discarded, what represents a waste of computing resources.The proposed method addresses this problem, filling the whole spectrogram with useful diagnostic contents.It relies on frequency shifting the Gaussian analysing window, as in sub-section 2.2,(13).If the current signal is low-pass filtered with a cut-off frequency f b (using, for example, a frequency filter as in[42]), then its spectrogram, built with a Gaussian window (5), will be non-zero only in the frequency band of diagnostic interest [0 − f b ].But, if the Gaussian window is shifted to a frequency f k = f b in (10), then the spectrogram of filtered current signal will appear in the frequency band [ f b − 2 f b ] (and will be zero outside this band).If the α parameters of both Gaussian windows (the shifted and the non-shifted one) are equal, the same spectral information will appear in both frequency bands.But, if the shifted Gaussian window has a different value of α, then the spectrograms in the frequency bands [0 − f b ] and [ f b − 2 f b ] will be different, each one corresponding to a different value of the parameter α.

Figure 1 .Figure 2 .
Figure 1.Multi-band window designed in Section 3 in the time domain (real part, top, imaginary part, middle), and in the frequency domain (bottom).This windows contains 20 Gaussian windows, located in 20 adjacent frequency bands of 100 Hz, spanning linearly the range [α = 1α = 700].

Figure 3 .
Figure 3. Synthetic current signal (17) used for illustrating the application of the proposed method.

Figure 4 .Figure 5 .
Figure 4. Multi-band window(15) in the time domain(top), and in the frequency domain (bottom).This windows contains 2 Gaussian windows, in 2 adjacent frequency bands of 100 Hz, one with α = 1 and the other one with α = 12.6.2.The spectrogram of the current signal (17) is built, using the multi-band window(18) as sliding window (Fig.5).It shows two elementary spectrograms in adjacent TF regions, obtained with a single run of the STFT algorithm.The bottom one (α = 1) locates the sinusoidal component at 50 Hz, but blurs the chirp component.On the contrary, the top one (α = 12.6) locates the chirp component, but widens the sinusoidal component.This second Gaussian window, and the spectrogram that it generates, have been shifted to the frequency band [100 Hz-200 Hz].

3 .
All the stacked, elementary spectrograms obtained in step 2 are shifted back to the frequency band [0-f b ], as shown in Fig.6.This process has a negligible computational cost, just the renumbering of the frequency axis of each elementary spectrogram.

Figure 6 .
Figure 6.Relocation in the frequency axis of the elementary spectrograms, so that all of them span the same frequency band [0-f b ].

4 .
The points with the same time-frequency coordinates in all the relocated spectrograms obtained in step 3 (Fig.6, are combined to give a single high resolution spectrogram of the TF region of diagnostic interest, in the frequency band [0 -f b ] (Fig.7).The combination process used in this work consists in selecting, for each point of this region, the minimum value obtained among all the relocated spectrograms.The final result shows with a high resolution both the sinusoidal component at 50 Hz and the chirp component.

Figure 7 .
Figure 7. High resolution spectrogram of the current signal (17): for each point of the TF region of interest, the minimum value obtained among all the stacked spectrograms of Fig. 5, right, is selected.The final result shows with a high resolution both the sinusoidal component at 50 Hz and the chirp component.

Figure 8 .
Figure 8. Start-up current of the motor of Appendix A with a broken bar fault.

Figure 9 . 2 .Figure 10 .
Figure 9. Multi-band window (21) in the time domain (real part, top, imaginary part, middle), and in the frequency domain (bottom).This windows contains 40 Gaussian windows, located in 40 adjacent frequency bands of 125 Hz, with values of α ranging from α min = 1 to α max = 700.2. The spectrogram of the current signal of Fig. 8 is built.First, the current signal is low-pass filtered, keeping only the frequency bins of its spectrum lower than 125 Hz.After, and using the multi-band window (21) as sliding window, the STFT algorithm (4) is applied, which generates the spectrogram sown in Fig. 10.This spectrogram contains 40 elementary spectrograms in adjacent TF regions (Fig. 10, right), obtained with 40 different Gaussian windows (Fig. 10, left), at the cost of a single run of the STFT algorithm (6 seconds with the computer of Appendix C).

Figure 11 . 4 .
Figure 11.Zoom of Fig. 10 showing two of the individual Gaussian windows contained in the multi-band window of Fig. 10 (left), and the corresponding spectrograms generated with them (right).The two zoomed bands corresponds to the spectrogram located in base frequency band [0-125 Hz] (bottom), which defines clearly the fundamental component, but blurs the fault harmonics, and to the spectrogram shifted to the frequency band [1500-1625 Hz] (top), which defines clearly the fault harmonics, but widens the fundamental component.

Figure 12 .
Figure 12. High resolution spectrogram of the current of Fig. 8: for each point of the TF region of interest, the minimum value obtained among all the stacked spectrograms of Fig. 10, right, is selected.The final result shows with a high resolution both the sinusoidal component at 50 Hz and the LSH fault component.

Table 1 .
Characteristic Frequencies of Different Types of IM Faults