Time Domain Analysis and Spectral Methods for Determining Rotational Speed of Rotary Machines

Recent Innovations in Mechatronics (RIiM) Vol. 3. (2016). No. 1-2. DOI: 10.17667/riim.2016.1-2/8. Time Domain Analysis and Spectral Methods for Determining Rotational Speed of Rotary Machines Gábor Gárdonyi Department of Mechatronics Optics and Engineering Informatics, Budapest University of Technology and Economics Budapest, Hungary [email protected]me.hu Krisztian SAMU PhD Department of Mechatronics Optics and Engineering Informatics, Budapest University of Technology and Economics Budapest, Hungary [email protected]me.hu Abstract — Accurate estimation of rotational speed of rotary machines has usually high priority in technical applications. This information should be calculated for many diagnostic algorithms, control or regulation processes. Incorrectly estimated values could occur serious disturbances in the operation of machines. Additional instrumentation often may be obstructed due to lack of space, but the construct of the machine may also affect the accuracy of measurement. In such cases, vibration diagnostic tools can be the disposal of difficulty. Mounting an acceleration sensor onto the outer surface of the measured device is not a major challenge. In most cases using time, frequency or quefrency domain analysis, it is possible to estimate the rotational speed of the analysed rotary machine. The calculated spectra and cepstra can help to determine the rotational speed more easily and more accurate than the time domain methods. This paper presents the comparison of these methods in terms of their usability and rotational speed estimation accuracy. A possible error of traditional optical measurement due to misalignment and benefits of the other methods are illustrated in this article via measured data series of a Brushless DC (BLDC) motor driven system. Keywords — cepstrum, diagnostic, rotary machine, rotational speed, spectrum I. INTRODUCTION This paper is concerned with the development of methods to calculate rotational speed of rotary machines using encoder or vibration signals. Rotary machines are an integral part of our everyday life, even if we do not perceive it directly in all cases. But if we think about it: rotating machines are in our vehicles, in many household appliances, children’s toys, modern robots, production equipment, and even mobile phones contain rotating parts. Rotational speed should ideally be measured for a high number of diagnostic algorithms, control or regulation processes. Incorrect values could cause serious disturbances in the operation of machines. There are plenty of common rotational speed measuring solutions to solve this problem. One can easily use electrical, optical principle based techniques, robust magnetic [1] or combined [2] techniques as well. Other specific solutions which have been published to date include image processing techniques [3] and electrostatic sensors [4]. The question may arise how one could measure the key parameters of an operating machine, particularly the rotating speed, if the engineers have not performed these measurements in the early stages of design. Mounting sensors into a device at a later stage is usually not an easy task, especially if it requires the modification of a rotating component – gear transmission system or rotor – which cannot be observed directly by the user. The difficulty is mainly caused by the compactness of modern devices. These parts are usually not directly accessible to users and there tends to be insufficient free space near the drive. In these instances, vibration diagnostics that uses additional acceleration sensor can be the best alternative solution. Vibration signal analysis is nowadays the most commonly used method for the troubleshooting and the condition monitoring of rotary machines. In recent years, research on vibration based rotational speed determination methods has become very popular. Many applications, further researches and alternative calculation methods can be found in the following references [5]-[10]. II. VIBRATION MEASUREMENT OF ROTARY MACHINES Mounting accelerometers onto the device surface takes not a major challenge. Furthermore there are also available non-contact vibration velocimeters, which have much less installation requirements. Nowadays in the case of most vibration measurement solutions piezoelectric accelerometers are used because of their small physical size, the good sensitivity and their especially wide dynamic range. Probably the main reason of this choice could be that accelerometers do not require an independent reference point. Instead of fixing it onto the surface, these sensors can also be used via special probe-tip. If more accurate measurement is required, magnetic fastening or screwing is recommended. If the surface is unsuitable for these type of fixings, mounting by special beeswax can be also a feasible solution to the problem. Further information about accelerometers and their mounting can be found in [11]. The rotational speed can be easily obtained from vibration signals that has been sampled in time. In addition it is also possible to get detailed information about condition of the measured system. Using advanced diagnostic the date of rotary machine breakdown can be well predicted as well. We can sort diagnostic algorithms into several categories: time domain analysis, frequency or quefrency domain analysis, Recent Innovations in Mechatronics (RIiM) Vol. 3. (2016). No. 1-2. DOI: 10.17667/riim.2016.1-2/8. higher order spectral analysis, order analysis, wavelet analysis and many kind of joint domain analysis methods as timespectra or time-cepstra. In this article the accuracy of estimated rotational speed calculated via methods in time and frequency domain and the cepstra are analysed. A. The analyzed system configuration The test system consisted of a BLDC servomotor (Faulhaber 3557-K024-CS) connected to a hysteresis braking motor (HB-20M-2). A wheel with radial oriented black and white stripes was fixed to the end of the assembly’s driven shaft. Thanks to this striped wheel the angular movement could be detected with a single reflective optical sensor. The piezoelectric accelerometer (PCB 356A33) was mounted onto the housing surface of the BLDC servomotor with special beeswax. III. POSSIBLE ERROR OF OPTICAL ROTATIONAL SPEED MEASUREMENT Considering the layout (Fig. 1), it becomes clear that mounting of an accelerometer takes much less time and effort as other alternatives. An advantage of the axle mounted striped wheel is that the rotational speed can be obtained with simplistic signal processing methods with good accuracy. Due to the mechanical structure – if the bearing is not suitable the system may leave the proposed operational range – intensive vibration of the analysed system could be observed. Because of this phenomenon, a rate of gross error of measurement has occurred. If the striped wheel have left the operational range of the used optical sensor, there would have been undetectable stripes. However, there are several correction methods which could be used to repair the calculated result. The result after the correction cannot be so accurate, as it would be if a better measurement assembly were used. Fig. 2 shows the fickleness of optical method above 2000 rpm (revolutions per minute) rotational speed due to the non-detected stripes of the wheel. In such instances the fault has caused a level of up to 60% relative error in the result. A section of the source signal from a steady state measurement in this high-speed range is interpreted in Fig. 3. This figure shows the measured acceleration and the voltage signal of the optical encoder in time. The Figure clearly shows that there were cyclically missing signal changes at a particular angular range. In these cases we can force correction techniques, but speed estimation based on vibration is also possible. Hiba! A hivatkozási forrás nem található. Fig. 1. Optical rotational speed measurement (left) and vibration acceleration measurement (right) of the analysed system configuration Fig. 2. Results of optical measurement of varying rotational speed Fig. 3. Measured signals of optical sensor (yellow) and the 3D accelerometer (blue, red, green) at steady state of the system Recent Innovations in Mechatronics (RIiM) Vol. 3. (2016). No. 1-2. DOI: 10.17667/riim.2016.1-2/8. IV. ROTATIONAL SPEED ESTIMATION BY TIME SIGNAL PROCESSING A. Comparison method of filtered time signal A special digital chart can be obtained from the sampled time signal. The source of cyclical changing vibration and the orientation of accelerometer must be taken into account at this calculation method. The calculated result describes the rotational speed. For detection of events in time signal, signal processing was performed in two main steps. Firstly smoothing was applied with a simple unweighted MA (Moving Average) low pass filter (Eq. 1). Eq. 2 represents the recursive formula of the unweighted MA calculation. Secondly we compared the calculated data series to a specified threshold level (Eq. 3). In the following equations xi is the ith sample of the input data series, yi is the ith sample of the output data and M marks the applied kernel size. TH represents the threshold level used for comparison. It is important to note that in the case of this experiment MA filtering has been enough to reach the expected results. In many cases running RMS (Root Mean Square) level calculation or other statistical values as skewness, kurtosis, etc. could be more effective.         1 0 1M j jix M iy (1)           1 0 1 1M j jix M iyiy (2)            THixif THixif iy 0 1 (3) The calculated results for the second measurement channel of the analysed assembly (vertical acceleration) is shown in Fig. 4. In this way we get a result that is similar to the optical encoder output. The disadvantage of this method is its low resolution. Usually there is only one or two detectable threshold crossing per revolution, so the speed can be determined just once or a few times per revolution. In addition, this method does not produce any information about the possible fluctuation within one revolution. In Fig. 4 every second rising edge, which occurs threshold crossing indicates a full period of operation. One full revolution takes approximately 13 to 14 ms time that equals to about 4500 rpm. B. Using Autocorrelation function for determining rotational speed The ACF (Auto Correlation Function) shows the averaged repetition time of the analysed section of the data series. The source of the repeated phenomenon can be an unbalanced part of the assembly, cyclically increased dissipative effect – such as bearing friction – or the uneven torque of the driving motor. All of these reasons produce periodic deviation in the measurable vibration acceleration. The correlation is a mathematical tool for finding repeating patterns, such as presence of periodic signal components obscured by random noise. ACF is the correlation of a signal with itself at different points in time. It shows how similar is the analysed data series to itself and gives information about period time of similarity. Prominent feature of ACF is that the calculation for periodic signal results periodic output. By contrast the uncorrelated components such noises have been eliminated by the transformation. The most commonly used estimate of the theoretical ACF of a WSS (wide-sense stationary) random process is the so-called biased estimate. It can be computed by the following equation:          1 0 1hN k XX hnxnx N hR (4) where   hRXX is the hth sample of the autocorrelation function, h is the lag and * denotes the complex conjugate. The input data is a length N realization of the random process. The main properties of the result of Eq. 4 are the followings:  The autocorrelation function is an even function.  The calculated value for   0 XX R is proportional to the signal energy.  Global maximum of the function is always in 0. Fig. 4. Smoothed acceleration signal (above) and the compared beat diagram (below) Recent Innovations in Mechatronics (RIiM) Vol. 3. (2016). No. 1-2. DOI: 10.17667/riim.2016.1-2/8. Fig. 5 shows the ACF that was calculated for the vertical vibration acceleration of the analysed system. The second prominent peak of the function after the zero location indicates the searched time period. An important property of this method is that it produces the averaged rotational speed for the related time section. This property causes that the method is effective just in case of the steady states of the examined system. A prominent advantage against the previously presented time based method is that the result is much less influenced by transient deviations and random noises. V. ROTATIONAL SPEED ESTIMATION IN FREQUENCY DOMAIN The traditional spectral analysis is concerned with the study of how the power of a signal is distributed in the frequency domain. After decomposing the original data series into sinusoidal components, it is relatively easy to detect a power content that corresponds to the rotational frequency [10]. The value of this peak on the X-axis represents the main operational frequency. In order to appreciate all of the detected frequency peaks one has to gain a thorough understanding how the analysed assembly works. In the case of a simple rotary machine the rotating speed can be determined directly from spectral results. Before any kind of frequency analysis is undertaken it is necessary to clarify which part of the assembly is the dominant source of vibration. The equation that describes rotating speed in rpm (Revolutions per Minute) is as follows: 1 60 fn  (5) where n is the rotating speed and f1 shows the main frequency in Hz. For example for internal combustion engines the crankshaft speed can be worked out according to the first main harmonic order when the number of cylinders is known. Assuming fi/2 is the frequency of the first main harmonic order k=i/2, f1 is the rotating frequency of the engine, the engine speed n will be [8]: ifkfn ii /120/60 2/2/  (6) A. Spectrum analysis Maybe the most commonly used method in this category is the so-called Auto Power Spectrum (APS) calculation. To compute the APS, the DFT (Discrete Fourier Transformation) of the signal is computed, and then multiplied by its complex conjugate [15]. Hence the magnitude of an APS is equal to the square magnitude of a DFT. This is one methodology, but APS can be calculated at several ways. The APSD (Auto Power Spectral Density) is the APS normalized to a 1 Hz bandwidth. That means APSD is the APS divided by the interval between frequency data points. The APSD can be also produced as the DFT of the ACF. The equations of DFT and IDFT (Inverse Discrete Fourier Transformation) are the followings:              1 0 2 1N n kn N i enx N kXnxDFT  (7)             1 0 2 N k kn N i ekXnxkXIDFT  (8)                         1 0 1 0 2 sin 2 cos N n N n kn N nxikn N nxkX  (9) The formula for power spectrum calculation is as follows:           22 ImRe kXkXkP  (10) where   kP is the APS,   kX is the Fourier transformed   nx signal and f represents the frequency in Hz. During spectral analysis, firstly signals are usually cut into smaller pieces, where overlapping and windowing can be applied. After the spectral calculation has been completed for every separate data series, it becomes possible to draw more reliable conclusions based on the averaged result. Alternatively, spectral analysis can also be calculated without averaging. However, averaging is a very effective method for Fig. 5. ACF of the measured vibration acceleration Recent Innovations in Mechatronics (RIiM) Vol. 3. (2016). No. 1-2. DOI: 10.17667/riim.2016.1-2/8. eliminating random noises in case of stochastic processes and getting a smoother output. In the following study 8192 wide Hanning window has been used for computing the averaged spectra. For the resulting plot, see Fig. 6. As can be seen in Fig. 6, the method gives us a representative result. After the DC component (0 Hz) the first prominent spectrum peak belongs to the base harmonics of vibration. The location of this peak on the frequency axis is the rotational frequency. The spectral analysis gives a result with linear frequency resolution. In this case the measurement data series has got a sampling rate of 51.2 kHz, the spectrum has been calculated with 8192 samples long window size. These details result a frequency resolution of Δf=6.25 Hz. The location of the peak is at 75 Hz that means a rotational speed of 4500 rpm. However the resolution causes 8.3% relative uncertainty. The reliability becomes better by higher speed. Spectrum based estimation is possible based on upper harmonics as well. Accordingly the location of a specific harmonic component can be defined with better relative accuracy. From this frequency value the natural frequency can be calculated by a simple division. The advantage of the upper harmonic based method has been the much lower reading error, however this method has got its disadvantages as well. As higher the upper harmonic’s frequency as lower the SNR (Signal to Noise Ratio). It has resulted problematic harmonic peak detection. The measured noise type can consider as white noise. It means that every frequency component is equally contained regarding the energy of noise. In spite of this the energy level of vibration components from mechanical operation is as lower as higher frequency is examined. The reasonable consequence of the previously described facts is the increasingly low SNR. B. Analysis by Short-Time Furier Transformation As mentioned above, the output of spectrum and cepstrum can be smoothed using average calculation methods. Whilst linear weighting usually produces a better result during a steady state analysis, logarithmic weighting has been much preferred for unsteady states or for real-time applications. Using another method, the spectra or cepstra is calculated for every windowed signal section as opposed to averaging. This partial result can be used also for analysing the coefficients’ changes over time. This consideration has brought to the Short-Time Fourier Transformation (STFT). Eq. 11 demonstrates the calculating method of STFT for discrete time signals.                m km N j emnwmxknXnxSTFT  2 , (11) where   nx is the source signal,   mnw  represents a real, even window and   knX , is considered to be the Short-Time Fourier Transform. The disadvantage of the above method is that, lacking averaging, the result is noisier than the traditional spectral calculation’s output. Nevertheless, the great advantage of STFT is that it gives evaluable results even in the case of variable rotational speed. In such instances the averaged spectrum produces smeared peaks. These can be more distorted because of the moving upper harmonics and other noises. The interpretability of the results obtained through STFT is always dependent on a compromise. With a large Fig. 6. APS of the measured vibration acceleration signal Fig. 7. Result of Short-Time Fourier Transformation Recent Innovations in Mechatronics (RIiM) Vol. 3. (2016). No. 1-2. DOI: 10.17667/riim.2016.1-2/8. window size we get a high frequency resolution. The larger the frequency resolution used, the more accurate the rotational speed which can be determined using the time-frequency representation. On the other hand, when using a small window size the resolution along the frequency axis becomes lower. However, the localizability over time becomes better. It becomes possible to observe or analyse transient state processes due to the better resolution over time. More detailed information about STFT and the results of further improvement researches have been published under the following references: [12], [16]-[18] The output of STFT is a two-variable distribution in the time-frequency domain. We can display this two-variable data series as an image and also post-process this result with common image processing methods. Thus we can detect ridges, edges or apply digital convolution filters. The output of ridge detection is shown in Fig. 8. The first ridge is at 75 Hz, which equals to 4500 rpm rotational speed. The applied window length and the outputs resolution is the same as in the case of the further presented averaged PSD method. C. Rotational speed estimation via Cepstrum analysis Cepstrum was originally defined as the power spectrum of the power spectrum’s logarithm. Later, a newer definition was coined; cepstrum being the inverse transform of the power spectrum’s logarithm [13]-Hiba! A hivatkozási forrás nem található., as expressed mathematically as:       kPIDFTnc log (12) In the literature four basic kinds of cepstral representations can be found. These are the real-, complex-, powerand phase cepstrums.  Real cepstrum:   ][log][ kXIDFTnc  (13)  Complex cepstrum:   ][log][ ˆkXIDFTnc  (14)  Power cepstrum:   2 2][log][ kXIDFTncp (15)  Phase cepstrum:                  ][Re ][Im arctan][ kX kX IDFTncF (16) A detailed study of the calculation methodology used for discrete data series and the particulars of the cepstral results can be found in [13]. For a number of further studies on the history and the application of cepstrum analysis could be found in Hiba! A hivatkozási forrás nem található.- [15][18]-[21]. Fig. 9 represents the real cepstrum of the analysed data series. The illustration shows the outputs calculated for all 3 channels of the 3D accelerometer. This method converts the signal into so-called quefrency components. One of the earliest applications of the cepstrum theory can be found in the study of signals containing echoes Fig. 8. STFT representation after image processing Fig. 9. Result of Cepstrum analysis Recent Innovations in Mechatronics (RIiM) Vol. 3. (2016). No. 1-2. DOI: 10.17667/riim.2016.1-2/8. and in the study of speech analysis. The application discussed in these studies was aimed at detecting the harmonic structure of measured sound. Such harmonic structures so-called harmonic families can be detected by gearbox analysis or alternatively by any kind of analysis of rotary machines. The cepstrum calculation of a signal only results an evaluable output if its spectrum contains harmonic families, as it searches for periodic changes in the spectrum. For example, a clear one-component sinusoidal signal can be examined well by power spectrum; however, its cepstrum produces much less useful information. As it has been shown previously by spectral analysis, windowing and averaging methods should be used also by cepstrum calculation. The presented result (Fig. 8) was made with the same basic parameters as by the power spectrum calculation. It means 8196 sample was the width of the used Hanning window. In this case the first dominant peak after the noisy section close to 0 quefrency belongs to the rotational speed or the so-called main harmonic component. The quefrency location of this peak shows the period length of one full revolution in seconds. The peak can be found at 0.0133s along the quefrency axis and the data series has the same resolution in quefrency as the source signal had in time. It means a quefrency step size of 0.0195 ms because of the applied 51.2 kHz sampling rate. These properties result a possible frequency estimation error raised from the resolution of the discrete series. The step size between two data points can be expressed relative to the actual quefrency value along the x axis. In the case of the represented results the method gives the expected 4500 rpm rotational speed. At this point the relative step size of the series was 0.146%. By lower quefrencies the result is usually much noisier. That is why this rotational speed estimating method has become unreliable at very high speeds. D. Rotational speed estimation by Time-Cepstrum Analysis The relationship between APS and STFT has been discussed in the previous chapter above. By analogy, we can define an algorithm to examine the cepstrum changes over time. During this method, sections of the time signal are windowed and cepstrum is then computed for these pieces. If we use these partial results without averaging them, it is possible to represent the cepstrum changing over time in three dimensional graphs. These graphs are usually called cepstrograms. A great advantage of cepstrum or cepstrogram representation is that the resolution along the quefrency axis is not influenced by the window size. It solely depends on the sampling rate of the source signal. However, specifying the suitable block size is not a clearly definable task. If we choose an overly short section the period time falls outside of the result’s quefrency range. In this case evaluating the TCA (Time-Cepstrum Analysis) produced output is problematic. Increasing the applied window size has a negative effect on the localizability over time. The same mechanisms of action has been discussed for STFT above. The TCA method gives the same result by rotating speed estimation as the averaged cepstrum output. It means we get the expected 4500 rpm rotating speed in the examined steady state section with prominent relative accuracy. VI. CONCLUSION The presented methods for rotational speed estimation pose a great convenience if the examined rotary machine hasn’t got any built-in sensor for rotating speed measuring purpose. It is possible to complete the assembly with other sensors, but it causes often difficulty because of the complexity of devices. Furthermore the gross error of a widely used solution with optical sensor was tested as well. For such a problem many of vibration diagnostic methods can supply solution. The easiest way is to determine the periodicity of vibration signal directly in the time domain. The usability of signal comparison is better after signal conditioning. We can also determine periodicity from the calculated autocorrelation function. There are several methods to transfer signal into other domains using spectra or cepstra calculation. The relative step size of the discrete series by rotational speed estimation depends on the method, the applied block size, the sampling rate of the source signal and the actual rotational speed. All of the presented methods are able to show us the main properties of the operation. The spectral analysis has the better relative resolution as higher the main frequency is. However the cepstrum has an inverse behaviour. In case of both methods the SNR became worse at high speeds. There Fig. 10. Result of Time-Cepstrum Analysis Recent Innovations in Mechatronics (RIiM) Vol. 3. (2016). No. 1-2. DOI: 10.17667/riim.2016.1-2/8. are also available methods for speed estimation by varying rotational speed, these are the discussed STFT and the TCA. These solutions have almost the same properties about relative accuracy as the spectrum or the cepstrum calculations. It is possible to develop such hybrid methods in digital signal processing, which can combine the above discussed ways of calculation. Such a hybrid calculation method should be based on the observation of signal components, SNR level, estimated main frequency etc. The relative accuracy of several methods can be estimated from the result of previous iterations and the hybrid logic should make the decision, which theory could give the best result. References [1] C. Giebeler, D.J. Adelerhof, A.E.T. Kuiper, J.B.A. van Zon, D. Oelgeschläger, G. Schulz, Robust GMR sensors for angle detection and rotation speed sensing, Sensors Actuators A: Physical 91 (2001) 16–20. [2] Y.S. Didosyan, H. Hauser, H. Wolfmayr, J. Nicolics, P. Fulmek, Magneto-optical rotational speed sensor, Sensors and Actuators A: Physical 106 (2003) 168–171. [3] X.C. Zhang, J.L. Chen, Z.T. Wang, N. Zhan, R.C. Wang, Digital image correlation using ring template and quadrilateral element for large rotation measurement, Optics and Lasers in Engineering 50 (2012) 922– 928. [4] Wang, L.J., Y.H. Hu, Y. Yan, X.C. Qian, Rotational speed measurement using electrostatic sensors with single or double electrodes, 2nd IET Renewable Power Generation Conference (RPG 2013), 2013, 1-12. [5] M.L. Adams, Rotating machinery vibration: From analysis to troubleshooting, CRC Press, New York, 2000 [6] Y. Yang, X.J. Dong, Z.K. Peng, W.M. Zhang, G. Meng, Vibration signal analysis using parameterized time–frequency method for features extraction of varying-speed rotary machinery, Journal of Sound and Vibration 335 (2015) 350-366. [7] Z. Feng, F. Chu, M.J. Zuo, Time–frequency analysis of time-varying modulated signals based on improved energy separation by iterative generalized demodulation, Journal of Sound and Vibration 330 (2011) 1225–1243. [8] H. Lin, K. Ding, A new method for measuring engine rotational speed based on the vibration and discrete spectrum correction technique, Measurement 46 (2013) 2056–2064. [9] K. Rodopoulos, C. Yiakopoulos, I. Antoniadis, A parametric approach for the estimation of the instantaneous speed of rotating machinery, Mechanical Systems and Signal Processing 44 (2014) 31-46. [10] K.C. Gryllias, I.A. Antoniadis, Estimation of the instantaneous rotation speed using complex shifted Morlet wavelets, Mechanical Systems and Signal Processing 38 (2013) 78–95 [11] M. Serridge, T.R. Licht, Piezoelectric Accelerometers and Vibration Preamplifiers, Brüel & Kjaer (1987) [12] B. Boashash, Time-Frequency Signal Analysis and Processing – A Comprehensive Reference, Elsevier (2003) [13] R.B. Randall, Cepstrum analysis and gearbox fault diagnosis, Tech. Rep. 13-150, Brüel & Kjær Application Notes, Copenhagen, Denmark, 1981 [14] R.B. Randall, A history of cepstrum analysis and its application to mechanical problems, Surveillance 7, Chartres, France, October 29-30, 2013, 11-16. [15] B. Liang, S.D. Iwnicki, Y. Zhao, Application of power spectrum, cepstrum, higher order spectrum and neural network analyses for induction motor fault diagnosis, Mechanical Systems and Signal Processing 39 (2013) 342–360. [16] H.K. Kwok, D.L. Jones, Improved instantaneous frequency estimation using an adaptive short-time Fourier transform, IEEE Trans. Signal Processing 48 (2000) 2964–2972. [17] G. Meltzer, Y.Y. Ivanov, Fault detection in gear drives with non-stationary rotational speed – Part I: The time-frequency approach, Mechanical Systems and Signal Processing 17 (2003) 1033–1047. [18] J. Urbanek, T. Barszcz, J. Antoni, A two-step procedure for estimation of instantaneous rotational speed with large fluctuations, Mechanical Systems and Signal Processing 38 (2013) 96–102. [19] M. Chamay, S.D. Oh, Y.J. Kim, Development of a diagnostic system using LPC/cepstrum analysis in machine vibration, Journal of Mechanical Science and Technology 27 (2013) 2629-2636. [20] I. Paraskevas, E. Chilton, M. Rangoussi, The hartley phase cepstrum as a tool for signal analysis, ITRW on Nonlinear Speech Processing (NOLISP 07) (2007) [21] W. B. Collis, P. R. White, J. K. Hammond, Higher-order spectra: The bispectrum and trispectrum, Mechanical Systems and Signal Processing (1998) 12(3), 375-394.