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

Kechik D. A.; Aslamov Yu. P.; Davydov I. G.

УДК 621

D. A. KECHIK, Y. P. ASLAMOV, I.G. DAVYDOV

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

Belarussian state university of informatics and radioelectronics

Problem of estimation of variated frequency of components ofpolyharmonic signals has been arose. Three-dimensional time-frequency representation ofsignals is usually used to resolve this problem. But simple and reliable method of instantaneous frequency tracking is needed. Frequency tracking method based on estimation of shifts ofpeaks 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 offrequency 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 offrequency track estimation, but analysis of time-frequency representation 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 vibrational diagnosing of rotary equipment was checked out using spectral interference method.

Keywords: instantaneous frequency, wavelets, spectrogram, correlation function, frequency tracking, frequency estimation, fundamental frequency

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 Fig. 1. Components of

a) 6)

Fig. 1. Vibrational signal observed under widely varying speed conditions: a) - spectrogram and estimated frequency track; b) - spectrum of initial and resampled signals

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-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 saN

tisfy the model: s(t) = ^ An(t)cos(nra(t)t + n^(t)).

n=1

Let fundamental frequency (FF, shaft rotation frequency of vibration signal or the first harmonic otherwise) ra(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 f1, harmonic number k and frequency shift s(t):

log*™ (fk (t)) = logbase (s(t) " k " f1) =

= lOg*ase (S(0) + lOg*ase (k) + ^gbae (f1) = (1)

= Slog (t) + log*ase (k) + log*ase (f X

where base is empirically selected base of the logarithm. Then distance between IF of any two consequent components is preserved the similar:

logbase (fk+1 (t)) - logbase (fk (t)) = ^gbase (k) . (2)

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:

CCFf, fh, t, 9) =ff\S(t, f) -S(0, f-0)df, (3, a)

shg (t) = argmax ,CCF (f, fh, t, 0), (3, b)

where 9 is lag of correlation function, f...fh 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 slog(t) of FF at each i-th time moment were estimated as frequencies according to lags of maximums of CCF(f, fh, t, 9) according to (3b). Then logarithmic track was translated to linear using formula: s(t) = basiss,og(t) dfLog, where dfLog 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:

S%(t) = (s(t) -1) -100%. (4)

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:

SNRlocal,times (i) = if0 _ ,

[ V <N(f )>2] - 2 (5)

= [f=fo (f ) ] „ <Nlow>2-((i-') - fo) , (5)

[V<Skf0)>2] [ V <S(kf0) >2] k=1 k=1 0

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), < Nlow > 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 iopt = argmaxiSNRiocal,times (i) , band"

width is defined as f0 ... iopt • f0. 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):

SNR (3) - SNR (2) =

= (E + E2 + E3)/(Nlw -2f0)-

-(E\ + E2)/(NL-f0) = (6)

= (E3 -E2 -E)/(N2owf

Ei = (A - Now)2.

Expression (6) is positive if and only if E3 > E2 + E1, 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 < S(k -f0)) >= Ak = {1,0.8,0.7},k = 1...3, global SNR (relation of power of signal and uniformly distributed in 0. FN noise, FN = 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 t13 = {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(t1, f) and its CCF with S(t2, f) and S(t3, f) have been presented at Figure 2d. Maximum of correlation function indicates shift of the components measured in samples of frequency. Estimated logarithmic track slog(t) and track expressed in per cents 5%(t) are depicted at Fig. 2e, 2f. Variation of instantaneous frequency has been captured by the method.

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: s(t) = ft) - f(0)]/ft0>100 %. Accuracy of frequency shift estimation was evaluated as root mean

squared error: RMSE (x, X) = [1 / «V (x - Xt )2 ](1/2),

where xi 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 3 a. 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: si = sin(2nf(1+kit) t), ki is slope of the i-th signal. RMSE of FF shift

з Р™' ^^мМВ^^ННИП^пНн^Н^Ну^ВВЦР^

а)

b)

450 500 550 600 650

Log Frequency, log(Hz)

c)

d)

e) f)

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

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

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

d)

c

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

c) d)

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

predicted increasing of error due to extension of bandwidth has been confirmed. Increasing of RMSE of shift has been shown for two values of frequency deviation speed (k = {1, 3}) at various SNRs.

Proposed method was also tested for signals with sinusoidal frequency modulation law (modulation frequency was 3 and 4 Hz) in the similar terms of accuracy. Dependence of RMSE of frequency shift was presented at Figure 3 c at different frequency and depth of frequency modulation. Increasing any of frequency and depth leads to significant increasing of error at various SNRs.

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 interpl 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 30f

(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.

REFERENCES

1. Barysenka, S. Y. Single-channel speech enhancement using inter-component phase relations / S. Y. Barysenka, V. I. Voro-biov, P. Mowlaee // Speech Commun.- 2018.- Vol. 99.- P. 144-160.

2. Barkov, A. V. Monitoring and diagnostics of rotary equipment relying on its vibration / A. V. Barkov, N.A. Barkova, A. Yu. Azovtsev.- Saint-Petersburg: publishing center of SMTU, 2000.- 159 p.

3. Kosmach, N. V. Approach of vibrational diagnosing of rolling bearing / Н. В. Космач, Ю. П. Асламов.- 2020.

4. Influence of changes in shaft rotational speed of rotary equipment on frequency-domain processing / Yu. P. Aslamov [et al.] // Doklady BGUIR.-2018.- Vol. 113, № 13.- P. 13-18.

5. Zhang, X. A new time synchronous average method for variable speed operating condition gearbox / X. Zhang, G. Wen, T. Wu // J. Vibroengineering.- 2012.- Vol. 14, № 4.- P. 1766-1774.

6. Aherwar, A. Vibration analysis techniques for gearbox diagnostic: A review / A. Aherwar, S. Khalid // Int. J. Adv. Eng. Technol.- 2012.- Vol. 3.- P. 4-12.

7. Bechhoefer, E. A Review of Time Synchronous Average Algorithms / E. Bechhoefer, M. Kingsley.- 2009.- P. 10.

8. Li, H. Order Bi-spectrum For Bearing Fault Monitoring and Diagnosis Under Run-up Condition / H. Li // J. Comput.-2011.- Vol. 6, № 9.

9. Birchmeier, J. R. Order tracking signal sampling process / J. R. Birchmeier. - 2002.

10. Rotating speed tracking and sampling method of variable speed mechanical fault diagnosis / Tang Deyao [et al.] - 2013.

11. Synchronous computed order tracking analytical approach of characteristic of rotating machines vibration signal / Wang Jiang [et al.].- 2015.

12. Jia Limin. The maximum power tracing controlling method that speed-changing oar-changing wind power generating set is followed the tracks of based on optimum resisting moment / Jia Limin, Liu Zhan, Lei Tao.- 2015.

13. Sparse wavelet decomposition of signals for solving vibration diagnostics problems / Y. Aslamov [et al.] // First World Congress on Condition Monitoring. - London, the British Institute of Non-Destructive Testing, 2017.- P. 11.

14. Kechik, D. A. Estimation of the instantaneous frequency of the spectral components of a nonstationary vibroacoustic signal by the segmented Prony method / D. A. Kechik // III All-russian acoustic conference / Politeh-press. - Saint-Petersburg, 2020.- P. 7.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

15. Kechik, D. Segmented Autoregression Pitch Estimation Method / D. Kechik, I. Davydov // 2020 International Conference on Dynamics and Vibroacoustics of Machines (DVM).- 2020.- P. 1-6.

16. Sparse Wavelet Decomposition with Redundant Dictionary for Vibration Waveform Analysis / Y. P. Aslamov [et al.] // 14th International Conference on Pattern Recognition and Information Processing.- Minsk, BSUIR, 2019.- P. 6.

17. Lin, S. Logarithmic Frequency Scaling and Consistent Frequency Coverage for the Selection of Auditory Filterbank Center Frequencies / S. Lin // ArXiv180100075 Cs Eess.- 2017.

18. Bendat, J. S. Engineering Applications of Correlation and Spectral Analysis / J. S. Bendat, A. G. Piersol.- Wiley, 1980.328 p.

19. Algorithms for refinement of the shaft rotational speed for solving the problems of vibration diagnostics of rotary equipment / Y. Aslamov [et al.] // First World Congress on Condition Monitoring.- London, the British Institute of NonDestructive Testing, 2017.- P. 11.

20. Time-Varying Autoregressions for Speaker Verification in Reverberant Conditions / V. Vestman [et al.] // Interspeech 2017: Interspeech 2017 / ISCA.- 2017.- P. 1512-1516.

ЛИТЕРАТУРА

1. Barysenka, S. Y. Single-channel speech enhancement using inter-component phase relations / S. Y. Barysenka, V. I. Voro-biov, P. Mowlaee // Speech Commun.- 2018.- Vol. 99.- P. 144-160.

2. Барков, А. В. Мониторинг и диагностика роторных машин по вибрации / А. В. Барков, Н. А. Баркова, А. Ю. Азов-цев.- СПб: Изд. центр СПбГМТУ, 2000.- 159 с.

3. Космач, Н. В. Способ вибрационной диагностики подшипников качения / Н. В. Космач, Ю. П. Асламов.- 2020.

4. Влияние изменения скорости вращения вала роторного оборудования на обработку в частотной области / Ю. П. Асламов [et al.] // Доклады БГУИР.- 2018.- Vol. 113, № 13.- P. 13-18.