=Paper= {{Paper |id=Vol-3039/short4 |storemode=property |title=Bi-periodically correlated random processes as a model for gear pair vibration |pdfUrl=https://ceur-ws.org/Vol-3039/short4.pdf |volume=Vol-3039 |authors=Roman Yuzefovych,Ihor Javorskyj,Oleh Lychak,George Trokhym,Pavlo Semenov |dblpUrl=https://dblp.org/rec/conf/ittap/YuzefovychJLTS21 }} ==Bi-periodically correlated random processes as a model for gear pair vibration== https://ceur-ws.org/Vol-3039/short4.pdf
                     Bi-periodically correlated random processes
                          as a model for gear pair vibration
Roman Yuzefovycha,b, Ihor Javorskyja,c, Oleh Lychaka, George Trokhyma, Pavlo Semenovd
a
  Karpenko Physico-mechanical institute of NAS of Ukraine, 5 Naukova Str., Lviv, 79060, Ukraine
b
  Lviv Polytechnic National University, 12 Bandera Str., Lviv, 79013, Ukraine
c
  Bydgoszcz University of Sciences and Technology, 7 Al. prof. S. Kaliskiego, Bydgoszcz, 85796, Poland
d
  Odessa National Maritime University, 34 Mechnikova Str.,Odessa, 65029, Ukraine
                Abstract
                The model of gear pair vibration in the form of bi-periodically correlated random processes
                (BPCRP) that describes its stochastic recurrence with two periods is proposed. Particular
                cases of this model are considered. It is shown that BPCRP model allows one to analyses
                unequally the mean and the covariance function of the additive and multiplicative
                components. There are considered technologies for the estimation of the Fourier coefficients
                of the mean and the covariance functions.
                Keywords 1
                Bi-periodically correlated random processes, gear pair, vibration.

1. Introduction
    The vibrations signals of the rotating elements can be characterized by their timely deviations
whose features are cyclic repetition and stochasticity. When faults arise in machinery, some non-
linear effects occur, and the interaction of different harmonics can be detected in vibration signal.
This interaction can be detected by the analysis of the parameters of the periodical (about periodical)
variation of the of the first and the second order moment functions of the random processes [1–4]
(they also are called periodically or about periodically correlated random processes [5–9]). Therefore,
it is preferable to select their parameters as the indicators for fault detection [10–17]. Gear pair excite
vibration signal because of two main reasons: the periodic deviation of teeth stiffness because of the
meshing phase and manufacturing inaccuracy. The manufacturing inaccuracy includes constant and
variable step errors of the teeth. The periodic deviation of the mesh stiffness results in the appearance
of the periodic components of the mesh frequency f m  rf 1  nf 2 and its multiples. Here f 1 and f 2 -
the rotation frequencies of the gear wheels and r and n are some natural numbers. The error of the
meshing step and the misalignment of axes and shafts are developed by the appearence of harmonics
with base frequencies equal to kf 1 and lf 2 and combination frequencies pf m  kf 1 , pf m  lf 2 , where
 p , k , l are an integer numbers. In addition, the direct spectra of vibration can include the
components that belong to some frequency band around the resonance frequency of the gear pair in
the case of a vibro-impact regime occurring.

2. Modeling of gear pair vibration
   The methods offered in [12, 18] for analysis of vibration of gear pair grounded on the transmission
error model considered in [19]:
                                                      
                        x   x e  W  x m   x 1   x 2   ,                          (1) 

ITTAP’2021: 1nd International Workshop on Information Technologies: Theoretical and Applied Problems, November 16–18, 2021,
Ternopil, Ukraine
EMAIL:       roman.yuzefovych@gmail.com;     javor@utp.edu.pl;    oleh.lychak2003@yahoo.com;      george.trokhym@gmail.com;
p.a.semenoff@gmail.com
ORCID: 0000-0001-5546-453X; 0000-0003-0243-6652; 0000-0001-5559-1969; 0000-0002-2472-1676; 0000-0003-4121-6011
           © 2021 Copyright for this paper by its authors.
           Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
           CEUR Workshop Proceedings (CEUR-WS.org)
where W is a some load and    t     is an angular position of the gear. The terms x   and     m

x   describe the contact properties of the gears, while terms x   and x   are caused by
 e                                                                                  1        2

manufacturing error. It is supposed that each term x   , i 1,2 is periodic with a rotation period
                                                                  i

Pi 1/ f i    of the corresponding gear. There are three periodic terms in (1), namely
x e   W  x m   , x e   x 1   and x e   x 2   , which are periodic functions with period
Pm 1/ f m , P1 and P2 . The model of the cyclostationarity offered in [12, 18] was obtained by
introducing the fluctuations of the angular position of the gears as a some random variable. The mean
function of this random process includes the harmonic with frequencies f m , f 1 and f 2 . The
covariance function consists of three different kinds of harmonics, in that, the harmonics with
frequencies that are a linear combination of the rotation frequencies kf 1  lf 2 , the harmonics of the
mesh frequency nf m and the harmonics with frequencies that are a linear combination of the mesh
frequency and the rotation frequencies, i.e. nf m  kf i . The first and the second order non-
stationarities have been substantiated by the processing of vibration signals measured on the gear
systems [12, 18], and the quantities that describe the structure of the cyclostationarity estimated by
means of synchronous averaging were proposed to be used for fault detection.
    In [20–22], after applying the synchronous averaging with the period P1 or P2 , the vibration signal
is expressed as

                                                                                   
                                   M
                          g t    Al 1  al t  cos 2 f l t  bl t   l ,                           (2)
                                   l 1

here M is the number of gear mesh harmonics, Al and l are the amplitude and the phase of the l th
harmonic respectively. The modulation effects are described by the functions 1  al t              and b t  ,
                                                                                                            l

which are periodic with the considered rotation period. These functions are closely approximate to the
signal’s deterministic component corresponding to one revolution of the selected gear.
   In [19, 23] the gear vibration signal is modeled as
                                                
                        x   x 1   x 2   x 1,2   x c   n  ,                             (3)
where x   and x   describe the deterministic periodic oscillations generated by the rotation of
          1           2

the output and input wheels respectively, x   is a component with period P  r P  r P ,
                                                        1,2                                      12   1 1       2 2

 x   is the second order cyclostationary process with period P and n   is a fluctuation
 c                                                                                      12
component. The deterministic part of the signal (3) can be extracted by means of synchronous
averaging with common period P12 of the shafts as far as it is possible [19]. Paper is dedicated to the
development of the cyclostationary models of gear vibrations considered in the literature, their
concretization, and the elaboration on this basis of other estimation techniques for the analysis of the
modulation effects occurring in the vibration signals as the faults originate.

3. BPCRP representation
   The efficiency of vibration signal processing for machinery condition monitoring can be explained
mostly by their possibility to develop modulations caused by the appearance of faults. The modulation
effects in the vibration model as a periodically correlated random processes (PCRP), which describe
the stochastic recurrence with one period can be explained by the jointly stationary random processes
     
k t in their harmonic representation [8, 9, 24]:
                                                                       2

                                           t    k t e
                                                                  ik        t
                                                                       P1
                                                                                ,
                                                 k Z
where Z is a set of integer numbers and P1 is a period of the rotations for one of the wheels.
Following this equation, we concludes that the modulation of the signals of two stochastic rhythms
provided by the rotation of two wheels can be explained as
                                                                                             2

                                                          t    k e
                                                                                        ik        t
                                                                               P2            P1
                                                                                                      ,                                    (4)
                                                                       k Z

where the harmonic of frequency 2 / P1 and its multiples are modulated for this once by PCRP with
                                               2

                       t    t e
                                          il        t
period P2 : k
              ( P2 )                           P2
                                   kl
                                                         .
                            l Z
   Then, for the random process (4), we have:
                                                         t    kl t e i kl t ,                                                     (5)
                                                                  k ,l Z

            are jointly stationary random processes and   k 2 / P   l 2 / P  . The random
where kl t                                                                                                 kl            1    2

processes presented by series (5) are bi-periodically correlated random processes (BPCRP) [9, 25,
26]. As we can see from (5), a BPCRP is a sum of the amplitude and phase modulated harmonics.
Here frequencies  kl are the linear combination of the two main frequencies 10  k 2 / P1 and                                      
01  l 2 / P2  . The modulating processes have the mathematical expectations mkl  E kl t 
which are the Fourier coefficients of the mean:
                                          m t   E  t    mkl e i kl t .                                                            (6)
                                                                              k ,l Z


                                                  
   For the covariance function R t ,  E  t  t   ,  t   t  m t , we have                                         
                                        R t ,   E  t    R ( )e ,                   kl
                                                                                                          i kl t
                                                                                                                                           (7)
                                                                              k ,l Z
where
                                               R kl     rp k ,q l , p ,q e
                                                                                                      i  pq 
                                                                                                                 ,                         (8)
                                                                  p ,q Z



                                                                      
and rpqkl   E  pq t kl t   ,  pq t   pq t  mpq are the cross-covariance functions of the
PCRP processes, and the “¯” signifies complex conjugation. Thus, the cross-covariance functions of
the modulating processes defines the Fourier coefficients of the covariance function (7) in which the
numbers are shifted by k and l . It follows from (8) those cross-correlations of modulating processes

kl t  with different numbers k and l lead to bi-periodical non-stationarity of the second order. As
the result of these correlations, it is appear the correlation of the spectral components, which can be
characterized by the appropriate Fourier transformation of (8):
                                                                         
                                           f kl                    Rkl  e i  d .
                                                                  1
                                                                     
                                                                 2 
                                                                                                                                           (9)

   It follows from (8) that
                                        f kl     f p k ,q l , p ,q    pq ,
                                                             p ,q Z
                                                                                                                
where
                                                                         
                                         f pqkl          rpqkl  e i  d ,
                                                        1
                                                           
                                                       2 
are the cross-spectral densities of the modulating processes  pq t . The functions (8) and (9) are                  
respectively the covariance and spectral components [9, 25, 26].
   The zeroth covariance component R00                                is determined by auto-covariance functions
rpq    E  pq t   pq t    : R00     rpq  e
                                                                             i  pq 
                                                                                         .
                                                  p ,q Z

   This covariance function of the stationary approximation for the BPCRP is averaged BPCRP
covariance function.
   The zeroth spectral component
                                                
                                 f 00    f pq    pq ,
                                                           p ,q Z
                                                                                       (10)   
is a power spectral density of the stationary approximation for the BPCRP. It defines the spectral
decomposition of the averaged in time instantaneous power R 0,t for the oscillations.                
   We should note that the covariance and the spectral components are the total characteristics of the
amplitude and the phase modulation of the BPCRP carrier harmonics. The zeroth spectral component,
as can be seen from (10), is a sum of the power spectral densities of the modulating processes  pq t                                       
shifted by  pq . The components f kl              explained in (9) are a sum of the shifted cross-spectral
densities for modulating processes. Their numbers differs by k and l . Proceeding from the above-
mentioned, it is possible to conclude that the zeroth spectral function f 00  describes the spectral                        
composition of the oscillations and the non-zeroth functions f kl  . It explains the correlations of the 
harmonics of this composition for the components with frequencies shifted by
                                       
kl  k 2 / P1  l 2 / P2 . When modulating processes of the corresponding numbers are
mutually correlated, than these correlations are not equal to zero.

4. Method for statistical analysis
    The time synchronous averaging (TSA) method was one of the early techniques used for the
analysis of hidden periodicities [27, 28]. If the hidden periodicity is presented and modeled as a
PCRP, then such technology was used for evaluation of its mean and covariance function [9, 25, 26].
It is so-called the coherent method [29, 30]. Synchronous averaging was also used for analysis of the
vibration signals, which are characterized by the recurrence of two or more periods [3, 7, 9, 13, 18,
22]. We consider below its application for the estimation of BPCRP characteristics.
    The coherent statistics of the BPCRP mean function have the form
                                                               N l 1

                                                                 t  nPl  ,
                                                           1
                                              mˆ l                                                                                        (11)
                                                       N l n 0
where Pl is one of the non-stationarity periods and N l is the number of averaged periods. The
mathematical expression of (11) for l  1 is equal to
                                                                              2
                                 N 1 1
                                                                                                                                P1 
             Emˆ 1 t           m1 t  nP1    m0k e
                                                                        ik          t
                             1
                                                                                           mkl e
                                                                               P2                              ik kl t
                                                                                                                          s N  l   ,
                            N 1 n 0                       k Z                              k Z l Z
                                                                                                                             1
                                                                                                                                P2 
where
                                                       
                                  P  i N 1 1 P2 lP1           P                P 
                            s N  l 1   e            sin  N 1 1  / N 1 sin   1  .
                                  P2                             P2               P2 
                               1




   If P1  nP2 and n is a natural number, then s N l P1 / P2
                                                                        1
                                                                                                   1 and Emˆ t   m t  , i.e. formula
                                                                                                                               1

(11) is the unbiased estimator of the BPCRP mean function. In other cases, formula (11) is a biased
estimator of the mean additive component with period P1 . The bias value depends on the ratio P1 / P2
and tends to zero as N 1  .
   Using (11), we can form the formulae
                                        P                         2                                P                         2
                                                            ik                                                         ik

                                           mˆ 1 t e                                                   mˆ 2 t e
                                                                       t                                                           t
                                  1 1                                                           1 2
                      mˆ k 0                                     P1
                                                                           dt , mˆ 0l                                        P2
                                                                                                                                       dt ,
                                 P1 0                                                           P2 0
which, in the general case, are the asymptotically unbiased estimators of the Fourier coefficients of
the mean additive components.
   It is easily see that unbiased estimators of the BPCRP mean function and its Fourier coefficients
can be obtained using synchronous averaging with common period P :
                                                                                                    P
                                              1 N 1
                           mˆ t                   t  nP  , mˆ kl                               mˆ t e i kl t dt .
                                                                                                1
                                                                                                P 0
                                                                                                                                                                  (12)
                                            N n 0
  Here N is the number of realization periods P which are averaged.
  Taking into account (11), we can form the coherent estimators of the covariance function and its
Fourier coefficients:
                      1 N
                                                              
          Rˆ t ,    t  nP  mˆ t  nP   t    nP  mˆ t    nP  ,
                      N n 0
                                                                                           (13)                                                           
                                                                           P
                                                   Rˆkl                    Rˆ t , e i kl t dt .
                                                                   1
                                                                   P 0
                                                                                                                                                                  (14)

   Using the synchronous averaging of the BPCRP samples over one of the periods P2 in the form
                            N2
         Rˆ t ,      t  nP2   mˆ t  nP2   t    nP2   mˆ t    nP2  ,
                      1
           l
                   N 2 n 0 
we can detect only the additive covariance components and determine their Fourier coefficients:
                                   P                              2                                         P                                2
                                                            ik                                                                         il
                 Rˆk 0             R1 t , e                        dt , Rˆ0l                      R2 t , e
                                                                       t                                                                           t
                                 1 1 ˆ                            P1                                    1 2 ˆ                                 P2
                                                                                                                                                       dt .
                                P1 0                                                                    P2 0
   Note that we must use in (13) the unbiased or asymptotically unbiased estimator of the mean
function.
   Component estimators are represented by trigonometric polynomials:
                                                                               L1
                                                       mˆ t    mˆ kl e i kl t ,                                                                              (15)
                                                                           k ,l  L1
                                                                             L2
                                                   Rˆ t ,    Rˆkl  e i kl t ,                                                                            (16)
                                                                       k ,l  L2

where Lr , r 1,2 are the numbers of the highest harmonics. The coefficients of the polynomials are
determined by the formulae
                                                                           T
                                                                                t e kl t dt ,
                                                                   1
                                                                  T T
                                                      mˆ kl                                                                                                      (17)

                                            T
                     Rˆkl            t   mˆ t   t     mˆ t    e
                                      1
                                   T T 
                                                                                                                              kl t
                                                                                                                                       dt ,                    (18)
                                          
where T is the length of signal realization. The number of harmonics to be taken into account in (15)
and (16) can be obtained on the basis of the results of experimental data processing by means of the
coherent method or stationary spectral estimation.
   In the general case, employing formulae (17) and (18) leads to an increase of the additional errors
caused by leakage effects. These effects are absent as T  NP . Formulae (17) and (18) can then be re-
written in the form of (11) and (13). Indeed,
                                       k  1 P
                          1 N 1
                                                                                            P
                                                                                                          1 N 1                                  
                                                    t e i kl t dt                                          t  nP  dt .
                                                                                    1
                mˆ kl        
                          NP k                                                      P0       e kl 
                                                                                               i  t

                                                                                                             N n
                                 0       kP                                                                    0                                
   Similarly,
                      P
                                         1 N 1                                                                                     
   Rˆkl                                     t  nP   mˆ t  nP   t    nP   mˆ t    nP  dt .
                 1
                       e kl 
                         i  t

             P0                0          N n                                                   
   The discrete estimators for the Fourier coefficients of the mean and covariance functions can be
formed by substituting the integral transformations (17) and (18) by integral sums:
                                              1 K 1
                                     mˆ kl    t e i kl nh ,
                                              K n 0

                                                                                                                      
                           K 1
           Rˆkl  rh      nh   mˆ  nh    n  r  h  mˆ n  r  h  e i kl nh .
                        1
                        K n 0                                                    
   Here P1  M 1  1 h , P2  M 2  1 h and T  Kh , where K  rN M 1  1  nN M 2  1 .
   To avoid the aliasing effects of the first and the second kinds [32], it is recommended to choose the
sampling interval h in accordance with the inequalities
                                                               Pi           Pi
                                                    h               , h         , i 1,2
                                                             2L1  1      4L2  1
    If these inequalities are satisfied, the expressions (15) and (16) can be considered as the
interpolation formulae for the estimators. We should note that in the case of T  NP the component
estimators coincide with the estimators determined by the least squares (LS) method [9, 31, 32].
However, using the LS method allows one to avoid the leakage errors in general case. These errors
can be significant in cases when the values of the rotation frequency and/or their combinations are
close. To construct the LS estimators for the mean and the covariance function we rewrite the series
(6) and (7) in the form
                                                                       M1

                                                                              
                                               m t   m0   mlc cos l t  mls sin l t ,
                                                                       l 1
                                                                                                                     
                                                                    M2
                                    R t ,   R 0     R lc   cos l t  R ls   sin l t 
                                                                   l 1
                                                                                          ,
 where ml  ml l 
                          1 2
                              2
                                                        
                                ml  imls , Rl    Rl l    Rlc    iRls   ,
                              1 c
                                                             1 2
                                                                      1
                                                                      2
                                                                 2
                                                                     2
                          m0  m00 , R0    R00   , l   l j    , l 1  1, L1 , l 2  1, L2 .
                                                                                              j 1    Pj
                                                                
and N 1  2L1 L1  1 , N 2  2L2 L2  1 . The LS estimators for the Fourier coefficients of mean and
covariance function are defined as the quantities which provide the minimum values of the quadratic
functions
                                                                                                                                 2
                                                                                                                       
                                                            
                                                                   T                   M1
     F1 mˆ 0 , mˆ 1c ,..., mˆ Mc , mˆ 1s ,..., mˆ Ms
                                        1                1
                                                                              
                                                                                                                            
                                                                    t   mˆ 0   mˆ lc cos l t  mˆ ls sin l t   dt , ,
                                                                                                                          
                                                                                                                                         (19)
                                                                   0                  l 1


    F2 Rˆ0   , Rˆ1c   ,..., RˆMc   , Rˆ1s   ,..., RˆMs    
                                               2                                     2           
                                                                                                   2
                                   T
                                               M2
                                                                                                    (20)
                       t ,   R0     R l   cos l t  R l   sin l t    dt ,
                                        ˆ            ˆ c                 ˆ s

                      0                     l 1                                          
                                                                                         
where  t ,   t  mˆ t   t    mˆ t    . They are the solutions of the system equations
which represent the necessary conditions for the existence of the minimum of functionals (19) and (20):
                              F1         F1            F1
                                    0,          0,           0 , r  1, M 1 ,                      (21)
                             mˆ 0        mˆ r
                                              c
                                                       mˆ rs
                             F2    F2      F2
                                 0       0        0 , r  1, M 2 .                             (22)
                             Bˆ    Bˆ c
                                             Bˆ s
                                0        r        r
   The lag-dependent vanishing of the covariance function is the enough sufficient condition of the
mean square consistency of the Fourier coefficients for the mean function. It also can be indicator of
the asymptotic unbiasness of the estimators of covariance component. This condition is also sufficient
for the consistency of mean square of the covariance component estimators for Gaussian BPCRP. For
similar purposes, the series procedures were introduced, the latest of them are self-adaptive noise
cancellation [33] and spectral method [34]. The best result are obtained using time synchronous
averaging, however it requires a separate operations including individual resampling in each
considered case.

5. Conclusions
   The advantage of the LS estimators is the absence of the leakage effect. The possible bias of the
LS estimators can be caused only by the previous inexact estimation of the mean function. When the
realization length increases, values of the component estimators and the variances for the LS are
quickly drawing together. So, the LS estimation can be rated as the preferable technique for statistical
processing of the PCRP experimental time series.

6. References
[1] W. A. Gardner, Introduction to Random Processes with Applications to Signals and Systems,
     Macmillan, New York, 1985.
[2] W. A. Gardner, Cyclostationarity in Communications and Signal Processing, IEEE Press, New
     York, 1994.
[3] A. Napolitano, Generalizations of Cyclostationary Signal Processing: Spectral Analysis and
     Applications. IEEE Press, 2012. doi:10.1002/9781118437926.
[4] A. Napolitano, Cyclostationary Processes and Time Series: Theory, Applications, and
     Generalizations. Elsevier, Academic Press, 2019. 10.1016/C2017-0-04240-4.
[5] E. G. Gladyshev, “Periodically and Almost-Periodically Correlated Random Processes with a
     Continuous Time Parameter,” Theory Prob. Appl., vol. 8, pp. 173–177, 1963.
[6] Y. Dragan, I. Javorskyj, Rhythmics of Sea Waving and Underwater Acoustic Signals. Naukova
     Dumka, Kyjiv, 1980 (in Russian).
[7] Y. Dragan, V. Rozhkov, I. Javorskyj, The Methods of Probabilistic Analysis of Oceanological
     Rhythmics. Gidrometeostat, Leningrad, 1987 (in Russian).
[8] H. L. Hurd, A. Miamee, Periodically Correlated Random Sequences: Spectral Theory and
     Practice. Wiley, New York, 2007. doi: 10.1002/9780470182833.
[9] I. Javorskyi, Mathematical models and analysis of stochastic oscillations. Lviv, Karpenko
     Physico-Mechanical Institute of NAS of Ukraine, 2013 (in Ukraine).
[10] V. Y. Mykhailyshyn, I. M. Javorskyi, Y. T. Vasylyna, O. P. Drabych and I. Yu. Isaev,
     “Probabilistic models and statistical methods for the analysis of vibrational signals in the
     problems of diagnostics of machines and structures,” Mater. Sci., vol. 33, pp. 655–672, 1997.
[11] A. C. McCormick, A. K. Nandi, “Cyclostationarity in rotating machine vibrations,” Mech. Syst.
     Signal Process., vol. 12, no. 2, pp. 225–242, 1998. doi: 10.1006/mssp.1997.0148.
[12] C. Capdessus, M. Sidahmed, J. L. Lacoume, “Cyclostationary processes: Application in gear
     fault early diagnostics,” Mech. Syst. Signal Process., vol. 14, no. 3, pp. 371–385, 2000. doi:
     10.1006/mssp.1999.1260.
[13] J. Antoni, F. Bonnardot, A. Raad, M. El Badaoui, “Cyclostationary modeling of rotating machine
     vibration signals,” Mech. Syst. Signal Process., vol. 18, pp. 1285–1314, 2004. doi:
     10.1016/S0888-3270(03)00088-8.
[14] J. Antoni, “Cyclostationarity by examples,” Mech. Syst. Signal Process., vol. 23, pp. 987–1036,
     2009. doi: 10.1016/j.ymssp.2008.10.010.
[15] R. B. Randall, J. Antoni, “Rolling element bearing diagnostics – A tutorial,” Mech. Syst. Signal
     Process., vol. 25, no. 2, pp. 485–520, 2011. doi: 10.1016/j.ymssp.2010.07.017.
[16] R. Zimroz, W. Bartelmus, “Gearbox condition estimation using cyclostationary properties of
     vibration signal,” Key Engineering Materials, vol. 413–414, pp. 471–478, 2009. doi:
     10.4028/www.scientific.net/KEM.413-414.471.
[17] I. Javorskyj, I. Kravets, I. Matsko, R. Yuzefovych, “Periodically correlated random processes:
     application in early diagnostics of mechanical systems,” Mech. Syst. Signal Process., vol. 83,
     pp. 406–438, 2017. doi: 10.1016/j.ymssp.2016.06.022.
[18] Z. K. Zhu, Z. H. Feng, F. R. Kong, “Cyclostationarity analysis for gearbox condition monitoring:
     Approaches and effectiveness,” Mech. Syst. Signal Process., vol. 19, no. 3, pp. 467–482, 2005.
     doi: 10.1016/j.ymssp.2004.02.007.
[19] W. D. Mark, “Analysis of the vibratory excitation of gear systems: basic theory,” J. Acoustical
     Soc. America, vol. 63, no. 5, pp. 1409–1430, 1978. doi: 10.1121/1.381876.
[20] P. D. McFadden, “Detecting fatigue cracks in gears by amplitude and phase demodulation of the
     meshing vibration,” J. Vib., Acoust., Stress,. and Reliab., vol. 108, no. 2, pp. 165–170, 1986. doi:
     10.1115/1.3269317.
[21] P. D. McFadden, “Examination of a technique for the early detection of failure in gears by signal
     processing of the time domain average of the meshing vibration,” Mech. Syst. Signal Process.,
     vol. 1, no. 2, pp. 173–183, 1987. doi:10.1016/0888-3270(87)90069-0.
[22] G. Dalpiaz, A. Rivola, R. Rubini, “Effectiveness and sensitivity of vibration processing
     techniques for local fault detection in gears,” Mech. Syst. Signal Process., vol. 14, no. 3,
     pp. 387–412, 2000. doi: 10.1006/mssp.1999.1294.
[23] J. Antoni, R. B. Randall, “Differential diagnosis of gear and bearing faults,” J. Vib. Acoust.,
     vol. 124, no. 2, pp. 165–171, 2002. doi: 10.1115/1.1456906.
[24] I. Javorskyj, V. Mykhailyshyn, “Probabilistic models and statistical analysis of stochastic
     oscillations,” Pattern Recogn. Image. Anal., vol. 6, no. 4, pp. 749–763, 1996.
[25] I. Javorskyj, R. Yuzefovych, I. Kravets, I. Matsko, “Methods of periodically correlated random
     processes and their generalizations,” in Cyclostationarity: Theory and Methods. Lecture Notes in
     Mechanical Engineering, F. Chaari, J. Leskow, A. Sanches-Ramires (Eds.). Springer, New York,
     pp. 73–93, 2014. doi: 10.1007/978-3-319-04187-2_6
[26] I. Javorskyj, O. Dzeryn, R. Yuzefovych, “Analysis of mean function discrete LSM-estimator for
     biperiodically nonstationary random signals,” Math. Model. Comput., vol. 6, no. 1, pp. 44–57,
     2019. doi: 10.23939/mmc2019.01.044.
[27] I. Javorskyj, J. Leśkow, I. Kravets, I. Isayev, E. Gajecka-Mirek, “Linear filtration methods for
     statistical analysis of periodically correlated random processes – Part II: Harmonic series
     representation,”       Signal      Process.,       vol. 91,     pp. 2506–2519,       2011.      doi:
     10.1016/j.sigpro.2011.04.031.
[28] I. Matsko, I. Javorskyj, R. Yuzefovych, Z. Zakrzewski, “Forced oscillations of cracked beam
     under the stochastic cyclic loading,” Mech. Syst. Signal Process., vol. 104, pp. 242–263, 2018.
     doi: 10.1016/j.ymssp.2017.08.021.
[29] I. Javorskyj, V. Mykhailyshyn, “Probabilistic models and investigation of hidden periodicities,”
     Appl. Math. Lett., vol. 9, no. 6, pp. 21–23, 1996. doi: 10.1016/0893-9659(96)00005-5.
[30] I. Javorskyj, D. Dehay, I. Kravets, “Component statistical analysis of second order hidden
     periodicities,” Digit. Signal Process., vol. 26, pp. 50–70, 2014. doi: 10.1016/j.dsp.2013.12.002.
[31] I. Javorskyj, R. Yuzefovych, I. Matsko, Z. Zakrzewski, J. Majewski, “Coherent covariance
     analysis of periodically correlated random processes for unknown non-stationarity period,” Digit.
     Signal Process., vol. 65, pp. 27–51, 2017. doi:10.1016/j.dsp.2017.02.013.
[32] I. Javorskyj, R. Yuzefovych, I. Matsko, I. Kravets, “The stochastic reccurence structure of
     geophysical phenomena,” in Cyclostationarity: Theory and Methods II. / Eds. F. Chaari,
     J. Leskow, A. Napolitano, R. Zimroz, A. Wylomanska, A. Dudek, New York : Springer
     International Publishing Switzerland, Applied Condition Monitoring, vol. 3, pp. 55–88, 2015.
[33] R.B. Randall, N. Sawalhi, M. Coats, “A comparison of methods for separation of deterministic
     and random signals,” Int. J. Condition Monitoring, vol. 1, no. 1, pp. 11–19, 2011.
[34] P. Borghesani, P. Pennacchi, R.B. Randall, N. Sawalhi, R. Ricci, “Application of cepstrum pre-
     whitening for the diagnosis of bearing faults under variable speed conditions,” Mech. Syst.
     Signal Process., vol. 36, no. 2, pp. 370–384, 2013.