METHOD OF ESTIMATION OF FREQUENCY VARIATION RELYING ON ESTIMATION OF SHIFT OF SPECTRAL PEAKS

Problem of estimation of variated frequency of components of polyharmonic signals has been arose. Three-dimension - al time-frequency representation of signals is usually used to resolve this problem. But simple and reliable method of instan - taneous frequency tracking is needed. Frequency tracking method based on estimation of shifts of peaks of spectrogram has been proposed in this paper. It is assumed that shift of spectral peaks of components of signal is proportional to variation of fundamental frequency. Logarithmic scaling of time-frequency representation is used to make spectral peaks equidistant. Temporal dependence of shift of spectral maximums is obtained using correlation of windowed spectrum at the first frame and spectrum of signal in the current window. Then obtained track is translated in linear scale. Proposed method does not estimate values of instantaneous frequency or central frequency of signal component but estimates its variation. Advantage of the method is that it can estimate frequency track even if range of frequency variation and its central value are known roughly or unknown at all. Multiple components do not interfere to estimate fundamental frequency variation. Reduction of bandwidth is recommended to increase accuracy of frequency track estimation, but analysis of time-frequency representa - tion containing a few components is also possible. Dependency of performance of analysis of synthetic signals using the method on various signal to noise ratios under different conditions was estimated. Applicability of the method for vibration - al diagnosing of rotary equipment was checked out using spectral interference method.


Introduction
Processing of non-stationary signals is complicated because the most convenient representation of speech [1], vibrational and vibration envelope signals [2,3] is polyharmonic.Frequencies of components of machinery vibration at low frequencies and components of envelope of vibration are dependent unambiguously on shaft rotation speed [2,3], and their instantaneous frequencies (IF) are proportional too [4].Frequency track is used to resample signal with equiangular step [5] and obtain spectrums that are suitable to processing by convenient methods.Result of such operation was demonstrated at

54
ОБРАБОТКА ИНФОРМАЦИИ И ПРИНЯТИЕ РЕШЕНИЙ signal whose frequency was widely variated were smeared, components of resampled signal are almost discrete in frequency domain.
Hardware tracking techniques are usually used to eliminate speed varying [6][7][8][9][10][11][12].Time-frequency signal processing methods may be used to estimate frequency track relying on signal.The most obvious approach is to filtrate signal at spectral areas containing the most of energy of signal (e. g. peaks of scalogram, smoothed power spectrum [13]) and compute IF of selected component, but it has been shown to be ineffective [14] and it is not appliable when components are intersected in frequency domain due to high (two times and more) frequency deviation.Alternative approach is segmented autoregression method [14,15].It is reliable at low SNR and rapid variations of frequency.Drawbacks of this method are runouts and not smooth track, that may be partially compensated by median smoothing (Fig. 1, a) [15].

Proposed method
In this paper frequency tracking method based on estimation of shift of spectral peaks is proposed.The signal is assumed to be polyharmonic and satisfy the model: Let fundamental frequency (FF, shaft rotation frequency of vibration signal or the first harmonic otherwise) ( ) t ω be variated in time.Then we cannot restore it relying on spectrum of signal (Fig. 1b).If signal is supposed to be stationary at short window, we can obtain spectrum of signal at short windows that are generally weighted and overlapped.Short time Fourier transform is commonly used procedure.Time-frequency representation (TFR) may be obtained by equivalent methods (wavelet transform, filter banks).STFT and DFT-modulated filter banks were trialed in this paper (Fig. 3a).TFR was obtained using STFT as modulus of complex valued spectrogram computed using 'spectrogram' function of MatLab package.The other approach assumes decomposition of signal using DFT-modulated filter bank (inverse FFT of windowed complex spectrum).Window function was Fourier spectrum of MORL wavelet [16] with central frequency equal to center of each frequency channel [14,15].Then signal at each frequency channel was substituted by its envelope (modulus of analytical signal).
TFR may be treated as instantaneous spectrum: its samples S(t, f) at fixed time moment t evaluates amplitudes of each component at this moment (while S(f) denotes overall spectrum).If instantaneous FF is changed at the next time moment, all spectral peaks will shift proportionally [4].Frequency axis of TFR is recommended to be logarithmic scaled to make spectral peaks be equidistant [17].Then, multiplicative coefficients of frequencies of components become additive.One can express IF of any harmonic as product of initial value of FF f 1 , harmonic number k and frequency shift s(t): where base is empirically selected base of the logarithm.Then distance between IF of any two consequent components is preserved the similar: Finding maximum of cross-correlation function (CCF) is conventionally used method to estimate shift between two signals of similar shape.Shift of frequencies of components of signal at each time moment is proposed to be estimated as lag of CCF of spectrums at the first time moment and at the current moment: where θ is lag of correlation function, f l …f h are frequency limits of TFR, hat over variable means its estimated value.
Positions of maximums of correlation functions indicate shift of all peaks of amplitude spectrum in log-frequency domain expressed in samples and then variation of fundamental frequency.Values of logarithmic track s log (t) of FF at each i-th time moment were estimated as frequencies according to lags of maximums of CCF(f l , f h , t, θ) according to (3b).Then logarithmic track was translated to linear using formula: , where df Log is uniform spacing interval in logarithmic scale.Obtained track denotes variation of instantaneous FF, times.It is appropriate to get frequency track estimation in per cents: ) Relation of power of periodical component of analyzed signal to wideband noise power (Signal to Noise Ratio, SNR) enhances due to low correlations of noise [18].To maximize SNR, preliminary selection of band limits is necessary.For vibrational signals a few harmonics of shaft rotation frequency and gear mesh frequency are recommended to use because they are the most prominent usually.Then, spectrogram is computed at frequencies defined by analysis of spectrograms of signals ensemble.Spectral representation of ideal polyharmonic signal and.
Additive white Gaussian noise is uniform spectrum with weighted deltas at multiple frequencies.Then, noise power may be estimated as integral of squared noisy samples in assigned band, signal power may be computed as sum of squared differences between amplitude of spectral peak and averaged noise level in the peaks vicinity.Thus, we can roughly define appropriate band limits through probable range of deviation of FF harmonics.Let us express local SNR through harmonics number i: where <x> means ensemble average of x, x means temporal average of x, N(f) is noise power density obtained as smoothed power spectrum of signal and noise S(f), low N < > is averaged noise level at low frequencies.Local SNR evaluates relation of power of signal and noise after bandpass filtration.Then, optimal number of considered harmonics is , ( ) For the sake of simplicity, we can rewrite expression of local SNR through amplitudes of the constant harmonics and uniform noise and compare SNR(2) and SNR(3): Expression ( 6) is positive if and only if E 3 > E 2 + E 1 , then i might be minimized to enhance SNR.It is appropriate to assign wide frequency range if there is lack of a prior knowledge about probable shaft frequency range.

Results
Performance of the method was tested on synthesized polyharmonic signals whose fundamental frequency was variated.Amplitudes of components were 0 ( )) {1,0.8,0.7}, 1 3 The proposed method does not estimate instantaneous or central frequency.Then, performance of the method was evaluated as accuracy of frequency track estimation (Fig. 3) and detection rate of harmonics of resampled signal (Fig. 4).
Known frequency track of synthesized signal was translated to frequency shift in per cents: Accuracy of frequency shift estimation was evaluated as root mean squared error: where  i x is estimation of measured variable at the i-th experiment, x is the true value of measured variable.
Performance of the method relied on TFR based on STFT and filter bank was compared at Figure 3a.Length of the signal was 10 seconds, frequency of FF variated according to linear law with 10 % per second growth speed (from 100 till 200 Hz) according to expression: s i = sin(2πf(1+k i t) t), k i is slope of the i-th signal.RMSE of FF shift of filter bank solution was significantly increased at lower SNRs due to side lobes of CCF of instantaneous spectra, but filter banks are recommended to be used to decompose short signals.
Performance of proposed method for estimation of shift of fundamental frequency of signals with different speed of frequency growth (10 %, 20 % 30 % per sec) has been shown at The Fig. 3b.Length of the signals was 1 sec, initial frequency was 100 Hz.RMSE does not depends strong on speed of frequency variation, especially at high SNRs.It is notable that RMSE is significantly greater if deviation of FF is greater while its speed is the same at low SNRs and approaches half of deviation, but at high SNRs it does not depend on frequency deviation.According to (6), performance of the method decreased if bandwidth was extended.RMSE of FF variation was compared for two cases: when analyzed band included two harmonics (bandwidth 90…270 Hz) and three harmonics (bandwidth 90…410 Hz) (Fig. 3d).Theoretically Detection rate of three harmonics (Fig. 4a,  4c) and detection rate of the first harmonic (Fig. 4b, 4d) were calculated.Estimated shift of frequency was used to resample the signal to compensate frequency variation using interpolation procedure (function interp1 of MatLab package, 'pchip' method) [15].Sample points were computed as phase of harmonic signal whose frequency was variated according to estimated track.Query points were calculated as phase of signal with the similar initial frequency that was constant.Instantaneous frequencies of components of resampled signal were constant and were equal to their initial values.Detection of harmonics was carried out using spectral interference method [19].This method takes into account accuracy of frequencies of multiple components, splitting of the components and presence of various spectral peaks in the same analysis frame and has been used for evaluation of performance of the frequency estimation methods [15].The highest spectral peak has been searched in spectral window of 10 % width of each component.Each spectral component is accepted to be detected if peaks of amplitude amounted 66 % of maximum of analysis frame were not detected in the same window and local SNR ≥ 3 dB.This condition restricts splitting (amplitude and number of prominent side components) and harmonic and wideband noise influence.Three harmonics shall be treated as detected if all of them were detected separately, and accuracy of central frequencies of multiple components better than 1 %.Accuracy was checked out using spectral interference method.Even small FF sweeping leads to leakage or splitting of the components, especially at high frequencies.Three harmonics of shaft rotation frequency are desirable for diagnosing using frequency domain methods [2], and it requires higher accuracy of shift estimation.Detection of spectral peak of FF indicated that deviation of estimated track from real one is not significant, and then FF of resampled signal is almost stationary.Detection of multiple harmonics requires higher accuracy.In some cases, three harmonics are not detected by spectral interference method due to leakage and offset of peaks.
In practice of processing of vibrational signals two-three harmonics of shaft rotation and gear mesh frequencies are usually present.Tracking of components of vibration envelope is also appropriate, especially if energy of useful components is low.Many harmonic components related with shaft rotation frequency are present in envelope of vibration of bearings and gearings [2] till 30•f 0 (the third harmonic of ball pass outer ring frequency of 6213 bearing) and higher.

Conclusions
The method of instantaneous frequency tracking for processing of poyharmonic signals has been proposed.Its effectiveness was tested on synthesized signals.High performance of the method has been shown for various values of parameters of tested signal and SNR.Error of shift estimation grows when SNR decreased lower than -4 dB under conditions of any frequency modulation law for signals of 1 sec length.Rate of single harmonic detection decreased at the same SNR values.Detection of three harmonics using spectral interference method requires higher SNR and starts to decrease approximately at 1 dB.Unlike segmented Prony method, the proposed one cannot estimate absolute FF and compare errors of frequency estimation.RMSE of frequency estimation of Prony method remained less 2 % of central frequency under the similar conditions and SNR > -5 dB.Error grows faster when SNR decreased under estimated threshold.Narrow band of analysis is recommended rather than wide band.It may be appropriate to use weighted least squares approximation to eliminate outliers and decrease error assuming some model of continuous frequency variation (e. g. [20]) and consider information carried by all components as well as envelope of vibration.

Fig. 1 .Fig. 1 .
Fig. 1.Vibrational signal observed under widely varying speed conditions: a) -spectrogram and estimated frequency track; b) -spectrum of initial and resampled signals global SNR (relation of power of signal and uniformly distributed in 0…F N noise, F N = 1 kHz is Nyquist frequency) variated in wide range.Spectrogram of synthesized signal with linear and logarithmic scaled frequencies has been shown at Fig. 2a), 2b) respectively.Examples of temporal slices of spectrogram at three time moments t 1…3 = {0, 1, 5} sec are presented at Figure 2c.It is equivalent of spectrum of temporal frame of signal or values of instantaneous spectrum at fixed time moments.It is obvious that spectral peaks are almost equidistant and are shifted along y-axis.ACF of S(t 1 , f) and its CCF with S(t 2 , f) and S(t 3 , f) have been presented at Figure 2d.Maximum of correlation function indicates shift of the components measured in samples of frequency.Estimated logarithmic track ŝ log (t) and track expressed in per cents ŝ % (t) are depicted at Fig. 2e, 2f.Variation of instantaneous frequency has been captured by the method.

Fig. 2 .
Fig. 2. Illustration of FF tracking using log-scaled TFR: a) Filter bank TFR; b) Log-scaled TFR; c) slices of log-scaled TFR; d) CCFs of spectral slices; e) log-scaled frequency track; f) estimation of FF track

Fig. 3 .
Fig. 3. Dependencies of RMSE of shift of FF on SNR: a) -comparison of two methods of TFR; b) -accuracy on conditions of different speed of frequency variation; c) -accuracy on sinusoidal frequency modulation conditions; d) -accuracy at different bandwidth

Fig. 4 -
Fig. 4 -Dependencies of detection rate of harmonics of resampled signal on SNR: a) -detection of three harmonics of signals with linear modulation; b) -detection of the first harmonic of signals with linear modulation; c) -detection of three harmonics of signals with sinusoidal modulation; d) -detection of the first harmonic of signals with sinusoidal modulation