Vibration Signal Extraction Based on FFT and Least Square Method

Fast Fourier Transform (FFT), widely used in spectrum analysis, is a powerful processing vibration signal to obtain the signal amplitude, frequency, and phase. However, the discrepancy between the FFT derived values and the real values could be introduced due to spectral leakage and spectral interference. This inconsistency is prohibited in some applications, for example, the extraction of vibration characteristics for power machinery and failure diagnostics. Therefore, the methodology to obtain the frequency domain's exact characteristics becomes one of the most concerning topics in vibration and signal processing. In this paper, a newly developed iterative method is presented in detail for high-accuracy characteristic extraction of multiple frequencies periodic vibration signals based on FFT and least square method. In the situation where the attenuation signal is superposed onto the periodic signal, the accurate characteristics of these two signals are also obtainable. Besides, simulation examples are provided, showing that the proposed method can be applied to the single-frequency signal, multiple frequency signals (including the signals of which the adjacent frequencies are close), and attenuation signal. The experimental results show that using the data processing method in this paper to extract the attenuation signals' characteristics, and the fitted attenuation curves are in good agreement with the actual attenuation signals. The method shown in this paper can be used for precisely extracting the characteristics of the periodic signal and attenuation signals in engineering.


I. INTRODUCTION
Rotating machinery, such as steam turbines, gas turbines, or large generator, has a significant role in the modern industry [1]. The method for monitoring and diagnosing rotating machinery, such as temperature, pressure, vibration monitoring, etc., can reduce the risk of economic loss and security problems caused by component failure [2]- [3]. Vibration monitoring is one of the most widely used methods [4]- [6]. Since the time domain signal is usually too complex to obtain effective information, some frequency-domain methods, such as Fourier transform, spectral analysis, wavelet transform, Hilbert transform, etc. are typically used. Fast Fourier transform (FFT) is the fundamental method for signal processing and is well accepted in engineering practice ever since Cooley-Tukey proposed [7]- [8]. However, the FFT processed amplitude, frequency, and phase may deviate from the frequency components' actual characteristics in the time domain given by spectral leakage due to non-integer-period truncation. For instance, as shown in Ref. [9]- [11], the maximum relative amplitude errors are 36.4% and 15.3% when using a rectangular window and Hanning window, respectively. The maximum error of frequency is half of the spectral resolution when using different windows, while the phase's maximum error can approach 90°. Besides, spectral interference can occur when two adjacent frequencies are too close, worsening the relative error. As a result, the data calculated by FFT without additional processing cannot be used directly in situations where high-precision amplitude, frequency, and phase data are needed. To name but a few, the experimental study of identifying the dynamic characteristic coefficients in sliding bearings [12]- [15], high precision rotor dynamic balancing tests [16]- [18], and the tests using impulse response method [19] to obtain the inherent characteristics of rotating rotors all require accurate data processing method.
In order to solve the mentioned spectral leakage problem, plenty of relevant works have been done by numerous researchers. The available methods for spectral correction are generally categorized into four groups: interpolation, FT continuous zoom, spectrum energy, and phase difference [20].
The interpolation method, initially proposed by Rife and Jain, uses two spectral lines to correct the spectrum. It is suitable for signals processed with the rectangular window [21]- [22]. In the following decades, many scholars improved the method, D. Agrez [23] proposes the weighted interpolated DFT algorithm in which three maximum spectral lines are used. These improvements extend the application areas of the interpolation method and enhance accuracy [24]- [26]. The influence of noise is also analyzed [27]- [28]. It is widely used due to easy implementation and high accuracy [29].
The spectral leakage can be minimized through FT continuous zoom, in which the spectrums are corrected by improving spectral resolution and reducing the spectrum analysis error. Zoom-FFT, a method developed from FT continuous zoom, can extract frequency range with high-frequency resolution without increasing the sampling number [30]- [31].
The spectrum energy method, firstly proposed by Carlo Offelli and Dario Petri, corrects the spectrum through the characters that the barycenters of the power spectrum of symmetric window function approach to the origin of coordinate [33]- [34].
The phase differences of the highest spectral line of two-time FFT can be used to restrain spectral leakage. This method has a good anti-noise performance but does not work well with the signals which are close to the adjacent frequencies [11], [20].
The methods mentioned above are accurate when processing a single component frequency signal (or two adjacent frequencies are far enough), but the accuracy is reduced significantly when spectral interference is serious. Luo [29] compares the accuracy of various interpolation methods with the Hanning window. When the normalized frequency approaches 1, negative frequency increases the spectral leakage and the frequency error. Luo and Xie [20] propose a phase difference method, in which the frequency error is 0.1 spectral bin when two adjacent frequencies are three spectral bins apart with a minimum 3-term window. For the method shown in Refs. [35]- [36], the maximum relative error of amplitude is 3.12%, and the maximum error of phase is 1.91°. The accuracy depends on the frequency correction method. This paper presents a high accuracy extraction method based on FFT and the least square method. Amplitude, frequency, and phase of the signal are obtained by searching for the best evaluation criteria of fitting accuracy. Accurate signal characteristics of the attenuation signal and the periodic signal with two close adjacent frequencies can be obtained using this method. The method is validated by the given calculation examples and percussion experiments.

A. LEAST SQUARE METHOD
The least square method is deduced from the criteria that minimizing the residual sum of squares [18]. The curves, Several sampling data of the observed signal is assumed as follows: Sampling moment T is an N×1 matrix and defined by 1 2 ( , , , ) T N T t t t   (2) Observed data of corresponding moment is an N×1 matrix and defined by 1 2 ( , , , ) T (2), the value of curvature function at the corresponding moment, which is an N×1 matrix, is given by Design matrix F, which is an N×m matrix and made by measurement scheme. It is given by 1 The coefficient matrix is defined by (6), which is an m×1 matrix, 1 2 ( , , , ) T m A a a a   (6) Based on (5) and (6), Ŷ can be expressed bŷ Y FA  (7) Residual is an N×1 matrix, which is given bŷ Y Y Y FA      (8) According to the criteria that minimize the residual sum of squares, the residual sum of squares is expressed as δ T δ and makes its partial derivative equal 0.
Based on (8) and (9), the coefficient matrix can be solved by Then curvature function can be obtained.

B. EVALUATION CRITERIA FOR CURVE FITTING
In general, there are three kinds of evaluation criteria for evaluating fitting precision. They are termed residual sum of squares, correlation coefficient, and standard deviation. The first two of them are adopted in this paper because the fitting effect of both the residual sum of squares and the standard deviation is almost identical. a: Residual sum of squares is defined by 2 2 ( ) i i y y     (11) A good curve fitting means a small δ 2 .
b: Correlation coefficient is defined by A good curve fitting means R is close to 1.
Where yi is the measured value; ŷi is the calculated value of the fitted curve; ӯi is the mean value of the measured value.

A. LIMITATION OF THE FOURIER TRANSFORM
The most typical signal in engineering is trigonometric polynomial. Its analytical expression is given by 2) The trigonometric formula (13.2) is used here in order to decompose the phase φ into the form of sine function plus cosine function. The design matrix can, therefore, be expressed concisely. The analytical expression of the signal is: Where m is the quantity of frequency component; fi is the frequency components of the signal; a2i-1 is the amplitude of the sine or cosine components.
An example of a single frequency signal is given by The time-domain figure of (15) is shown in Fig. 1. Another example of a double frequency signal is given by  Assuming the sampling frequency is 1024Hz, and the sampling number is 1024, (15) and (16) is processed using FFT. The frequency spectrum X(k) of each signal is shown in Fig. 3 and Fig. 4. There is no spectral leakage, and all the parameters of the analyzed signals, i.e., amplitude, frequency, and phase, are the same as the true signal. Assuming the sampling frequency is 1000Hz, and the sampling number is 1024, (15) and (16) is processed again using FFT. The frequency spectrum X(k) of each signal is shown in Fig. 5 and Fig. 6. Spectral leakage occurs, and the value of the spectral line is the envelope function value at the corresponding point. There are differences between the FFT results and the true signal. When the sampling time is not an integer multiple of the signal period, spectrum leakage will occur.
In addition, the spectral lines of the two frequency components interfere with each other, and the main peak became higher at a shifted location. This phenomenon is called spectral interference.
In general, if the total time is not the integral multiple of the period of the signal, spectral leakage will happen. Simultaneously, provided that two adjacent frequency components are too close, spectral interference will occur. It is hard to realize the above-mentioned conditions in signal processing. The character calculated by FFT is not accurate enough and will cause an error in the following calculation. Therefore, an improved method is needed to process the signal in order to obtain the right signal characteristics.

B. EXTRACTIONE METHOD FOR SIGNAL CHARACTERISTICS OF PERIODIC SIGNALS
In this section, we will introduce the method on how to use the precision evaluation criteria of the least square method to realize the least square fit.

1) SUM OF RESIDUAL SQUARES OF THE FITTING CURVE AND FREQUENCY ERROR
If there is only one frequency component (or the adjacent frequencies is far enough), the real frequency is located between the corresponding highest spectral line and the second-highest spectral line. If two adjacent frequencies are so close that the side lobes of two real frequencies interfere with each other, the real frequency may disappear from the frequency section bounded by the highest and the second-highest spectral line. In this case, an enlarged section is needed to retain the real frequency. In general, the real frequency is contained in the section given by Where fi is the i th real frequency component, and ki is the spectral line number.
The frequency converted from the highest spectral line is defined as the estimated frequency. It follows The difference between estimated frequency and real frequency is named frequency error, which is given by = The residual sum of squares of the curve is big, which is fitted by using the design matrix composed of the estimated frequency. The sampling signal is given by The design matrix is given by From (10) and (21), the residual sum of the squares of the curve is computed by The influence of the frequency error on the sum of squared residuals will be discussed next. Take the following signal as an example.
There are two frequency components in the signal, f1 = 10 Hz and f2 = 20 Hz. f1 is ranged and f2 is fixed. Using the least square method to fit the curve, the sum of squared residuals is calculated by using f2 and variable f1 (to simulate fe1). The relationship between frequency error and the sum of squared residuals is shown in Fig. 7. The difference between the real frequency f1 and the estimated frequency fe1 is the abscissa and the sum of squared residuals is the ordinate.   As is shown in Fig. 7, the residual sum of squares δ 2 = 0 in the circumstance of frequency error Δf = 0. With increasing frequency error, the sum of squares of residuals also increases. When Δf continues to increase, δ 2 will increase proportionally. When Δf increases to a certain value Δfm, δ 2 will reach the maximal value, and then fluctuation occurs. In section [0, Δfm], the relationship between Δf and δ 2 is a monotone increasing relation. At the same time, it is worth noticing that the value of Δfm changes with the variations of sampling frequency. Fig. 8 shows that Δfm and Fs have a linear relationship with different sampling number N. Fig. 9 shows that Δfm and N have an inversely proportional relation at different sampling frequency Fs. A conclusion that can be drawn from this is that / m s f F N   .

2) EXTRACTION METHOD FOR SIGNAL CHARACTERISTICS BASED ON THE LEAST SQUARE METHOD
Since the residual sum of squares, δ 2 will decrease with the diminution of Δfi in [0, Δfm], the real frequency fi can be approached by seeking the estimated frequency fei with the least residual sum of the square. The specific procedures are listed as follows: a: Applying FFT on the signal with sampling number N and sampling frequency Fs; spectral lines ki can be obtained to be the original estimated frequencies fei = ki×Fs/N. The estimated sections are listed as follows: Fs square fitting is done by using fe1 and another fei, and the residual sum of squares δ 2 in the corresponding different α is obtained. The smaller δ 2 is, the better α is. If α is good enough, the estimated frequency fe1 related to this α can approach real frequency f1. Replace the first estimated frequency by it. The method of searching the best coefficient α can be carried out by searching α related to the smallest and sub smallest δ 2 in the specified shrinking frequency interval. The interval is initially defined by   , ia ib f f and then is divided into several small intervals, 10, for instance. The new interval is determined by frequencies related to the smallest and sub smallest δ 2 in the last interval. c: Renew other frequencies in sequence, repeat step 2) and other real frequencies can be confirmed. d: After all frequencies are confirmed, repeat step 2)-3) to improve the accuracy of the result. Fig. 10 shows the calculation flow diagram. FIGURE 11. Attenuation signal added in a periodic signal.

FIGURE 12. Complex attenuation signal added in a periodic signal
The condition usually encountered is that the attenuation signal is added to a periodic signal when shock acts on a shaft in steady operation. The damping coefficient and vibration frequency can reflect the over damped oscillator and nature frequency well, so they are important parameters. It is assumed that the signal is measured from the moment T0, and the measured signal is periodic. After a period of time T, a shock occurs on the shaft, and a new signal component is added. The signal is an attenuation signal because of the damping coefficient α. The measurement stops at the moment T1.
During the time period [T0, T], the signal is expressed by During the time period [T, T1], the newly added signal is expressed by Where ε(t-T) is the step function, which is given by From (24) and (25), the signal is represented by

1) CALCULATION FOR COEFFICIENT MATRIX OF PERIODIC SIGNAL AND CONFIRMATION
Because (27) is superimposed by two functions, and the time of shock occurrence cannot be confirmed, the curve cannot be fitted directly. The function section is divided into two parts for discussion, i.e., [T0, T] and [T, T1]. In [T0, T], the method mentioned in section II. B. is used to fit the periodic signal.
In the general situation, Fig. 11, for example, the sampling number N1 belongs to [T0, T] can be obtained through observation. However, the sampling number of some complex signals, Fig. 12, for example, cannot be obtained through inspecting the length of the signal section. The method is introduced as follows.
Suppose a time section; there are Ni sampling points in the section. If the Ni points are in front of the moment T, the residual sum of square δ 2 calculated by the method mentioned in section II. B. will be close to 0 (the correlation coefficient R will approach 1). On the contrary, if some points out of the scope, δ 2 will deviate from 0 (R will deviate from 1).
Change Ni from the dimension of the coefficient matrix to a certain value. Using the method in section II. B. to calculate the coresponding δ 2 (or R), respectively. The derived curves are shown in Fig. 13 and Fig. 14.
As is shown in Fig. 13 and Fig. 14, before the N1 th point, the curve fits nicely. From the (N1+1) th point, the residual sum of square δ 2 and correlation coefficient R begin to deviate from the former value sharply. It means that the moment T is between the time represented by the N1 th point and (N1+1) th point. Since T is limited by [N1/Fs, (N1+1)/Fs], the T with higher accuracy can be obtained by improving the sampling frequency Fs. Then the signal in [T0, T] is confirmed by the method mentioned in section II. B.
After ye ʹ is divided by e -αtʹ , (28) is transformed as  It is a trigonometric polynomial on the right side of the equation and easy to be fitted. The only task left to do is to search for a right damping coefficient α to minimize the residual sum of square δ 2 . Assuming   c: If a higher high-precision α is needed, select α1 and α2 got from the last step to build a new section [α1, α2], repeat step 2)-3) to improve the accuracy.
As is shown in Fig. 15 and Fig. 16, the best fitting is at the extreme point of the δ 2 (or R) curve. Therefore, the value of α and the frequency is confirmed.

IV. CALCULATING EXAMPLES
In order to verify the correctness and accuracy of the signal feature extraction method introduced in Section II, this method is compared with the existing signal processing methods of rectangular window function and Hanning window function to highlight its advantages in signal feature extraction; some examples are designed in this section. These signals, which have 16 significant figures, are processed by FFT with the assumed parameters. The sampling frequency is 1000 Hz, and the sampling number is 1024. MATLAB is used for programming. Relative errors of frequencies and amplitudes are calculated through the program.

A. PERIODIC SIGNAL EXAMPLES
A single frequency signal is given by Iterations are set as aof = 20 and aof1 = 0; the results of the fitted frequencies of the signal frequency signal are listed in It is a quasi-periodic signal. In the signal, there are two adjacent frequencies, of which the frequency difference is 1.2 spectral bin. Iterations are set as aof = 20 and aof1 = 60; the results of the calculation are listed in In this signal, there are four adjacent frequencies, of which the difference is two spectral bins. Iterations are set as aof = 20 and aof1 = 5. The results are listed in TABLE  3.1 to TABLE 3.2.  TABLE 3.1 shows the results of the Fitted frequencies of  the multiple frequency signals; TABLE 3.2 shows the results of Fitted amplitudes of the multiple frequency signals. The results show that the relative error between the fitting value and the real value is smaller than that of the traditional rectangular window function and Hanning window function when the signal features are extracted by the proposed data processing method.
When multiple frequencies which are continuous and closed, fitting accuracy is lower than the results of the single frequency. The simulation shows that the residual sum of the square cannot be reduced due to the influence of the storage bit number. When one of the four frequencies is changed, the change in frequency causes a small change in the sum of squares of residuals. When the maximum number of storage bit is reached, the smaller residual variation cannot be reflected in the computer data, but from the relative error between the fitting results and the real results, the recognition accuracy can still meet the general engineering requirements.

B. ATTENUATION SIGNAL ADDED IN THE PERIODIC SIGNAL EXAMPLE
An attenuation signal added in the periodic signal is given by Iterations are set as aof = 20 and aof1 = 1, according to the method in section II. C. Points before 531 are periodic signals.
The results of the periodic signal are listed in

V. EXPERIMENTAL VERIFICATION
This section uses percussion experiments to verify the accuracy of the data processing method in this paper; the verification experiments are divided into rotor-bearing system percussion experiment and straight blade vibration percussion experiment.

A. VERIFICATION EXPERIMENT 1
Verification experiment 1 uses the rotor-bearing system experiment platform to carry on the percussion experiment to verify the accuracy of the data processing method in this paper. Fig. 17 shows the rotor-bearing system experimental platform. The experimental device includes the main body of the test bench, driving motor, signal acquisition, processing and analysis system. The main body and driving motor of the experimental platform are positioned by the slot slide rail and fixed on the cast iron platform by bolts. A single-wheel disc rotor is supported by two sliding bearings, and the sliding bearings are independently supplied with oil by oil cups. The diameter and thickness of the disc are 78 mm and 14 mm, respectively. The driving motor is a 22 Kw DC motor. The speed of the DC motor is controlled by the photoelectric speed sensor. The main body of the test bench and the driving motor are connected by elastic corrugated coupling. A set of eddy current displacement sensor is installed near the rotor disc to measure the vertical and horizontal displacement of the rotor. The displacement signals are displayed and recorded by the data acquisition system.
After the installation and debugging of the test bench, the percussion experiments are carried out on the rotor-bearing system test bench. When the rotating shaft is stationary, and the rotating speed is 500 rpm, the rotor disc is excited by the rubber hammer. The time-domain curves of the vertical displacement change of the rotor-bearing system after being excited are collected by the signal acquisition system. The sampling frequency is 16384 Hz, and the sampling number is 16384.
When processing the data signals collected by the rotor-bearing system test bench, Firstly, the vibration signal feature extraction method based on FFT and least square method described in this paper is used to identify the frequencies and amplitudes of periodic signals, By subtracting the identified periodic signals from the original time-domain signals, the time-domain signals with attenuation motions can be obtained. Furthermore, the attenuation coefficients and frequencies of attenuation motions are extracted, and the curve fitting is carried out. The fitting results are evaluated by coincidence between the fitting curves and the actual attenuation curves. Fig. 18 and Fig. 19 show the time domain curves of the rotor-bearing system in the vertical direction of the shaft after being excited by a hammer when the shaft is stationary, and the rotating speed is 500 rpm. The vibration signal feature extraction method based on FFT and least square method described in this paper is used to process the collected signals. Fig. 20 and Fig. 21 are the comparison diagrams of the actual attenuation curves and the fitted curves of the attenuation signals after the acquisition signals are processed. The red curves represent the fitting curves, and the black curves represent the actual attenuation curves of the attenuation motions. TABLE 6 shows the results of attenuation signals.
According to the time-domain curves of the signals collected from the two taps, the time-domain signals of the shaft before being excited by the rubber hammer is Periodic. When the shaft is excited by the rubber hammer, the vertical displacement increases and decays. Due to the large damping ratio of the system, the attenuation motion disappears quickly, and the time-domain signal of the rotating shaft returns to the periodic signal. Using the vibration signal feature extraction method based on FFT and least square method described in this paper to process the collected signals and intercept a section of the attenuation curves for feature recognition. TABLE 6 shows the rotating speed and attenuation frequencies of the identified attenuation signals; from Fig. 20 and Fig. 21, the attenuation curves are fitted by the attenuation coefficients of the identified attenuation signals, and the attenuation frequencies are in good agreement with the actual attenuation curves in the intercepted time periods.

B. VERIFICATION EXPERIMENT 2
Verification experiment 2 uses the blade vibration characteristics research experimental platform to perform percussion experiments to verify the accuracy of the data processing method in this paper. Fig. 22 shows the experimental platform for blade vibration characteristics research. The experimental equipment mainly includes the blade vibration excitation system, the blade vibration experimental platform, the signal acquisition and analysis system. The blade vibration test bench includes a straight blade and blade fixing devices. Straight blade with length, width, and thickness of 390 mm, 82 mm, and 2 mm are fixed on the blade vibration test bench. The electric vibration exciter provides the power source of blade vibration, and the exciting end of the vibration exciter is connected with the straight blade through a long screw. The acquisition probe of the accelerometer is adsorbed on the blade at a distance of 65 mm from the free end of the blade. The blade vibration signals collected by the acceleration sensor are displayed and recorded by the signal collection and analysis system. After the installation and debugging of the test bench, the percussion experiments are carried out on the blade vibration test bench. When the blade is stationary and the blade is forced to vibrate at 35 Hz, the blade is excited by the rubber hammer at the position 100 mm away from the free end of the blade. The time-domain curves of the acceleration change of the blade in the vertical direction after being excited are collected by the data acquisition system. The sampling frequency of the blade vibration test is 2048 Hz, and the sampling number is 6144.
Using the same data signal processing method as Section V. A, Fig. 23 and Fig. 24 show the time-domain curves when the blade is stationary, and the blade is forced to vibrate at 35 Hz. Fig. 25 and Fig. 26 are the comparison diagrams of the actual attenuation curves and the fitted curves of the attenuation signals after the acquisition signals are processed. The red curves represent the fitting curves, and the black curves represent the actual attenuation curves of the attenuation signals.  Fig. 23 and Fig. 24, after the blade is excited by the rubber hammer, it will decay motion under the action of system damping and restore the original motion state after the damping motion disappears. Using the vibration signal  feature extraction method based on FFT and least square method described in this paper to process the collected signals. When processing the attenuated signal, intercept the less clutter section of the attenuation curve for signal feature identification. TABLE 7 shows the vibration frequencies and attenuation frequencies of the attenuation signals identified after the blade is excited. From Fig. 25 and Fig. 26, the fitted attenuation curves and the attenuation curves of the actual attenuation signals are in good agreement in the intercepted time periods.

VI. CONCLUSION
A method for vibration signal extraction to acquire high-precision signal characteristics in the frequency domain is presented. Preliminary approximation information of vibration signal is obtained by Fast Fourier Transform (FFT) and Least Square Method is introduced to fit the curve and evaluation criteria. The residual sum of squares and correlation coefficient are used to evaluate the fitting results. The relationship between frequency error and residual sum of the square, which is used to evaluate the fitting results of the periodic signal, is discussed. Thus, an extraction method is proposed in which spectral correction of the periodic signal is carried out by searching for the best evaluation criteria. In the case of attenuation signal added in the periodic signal, the relationship between different estimated damping coefficients and the residual sum of the square of fitting results is also discussed. Then the method for obtaining the damping coefficient is proposed. The conclusions of this paper are as follows: Secondly, single, multiple-continuous, multiple, frequency signals, or periodic signals with added attenuated signals are analyzed through examples. High-accuracy signal characteristics are extracted, which proves that the method is correct and accurate. In addition, spectral interference is overcome while multiple continuous frequency signal is processed.
Thirdly, the experimental results show that the vibration signal feature extraction method based on FFT and least square method can solve the problems of spectrum leakage and spectrum interference caused by FFT, the attenuation curves fitted by the attenuation coefficients and attenuation frequencies of the attenuation signals identified by the method in this paper are in good agreement with the actual attenuation curves of the attenuation signals.