Robust coherent and incoherent statistics for detection of hidden periodicity in models with non‑Gaussian additive noise

,

scenarios.They are observed in mechanical systems [1,2], climatology, and meteorological processes [3,4].Furthermore, financial time series can also exhibit PC behavior [5,6].The concept of PC processes was initially introduced by L. Guzdenko and E. G. Gladyshev [7,8] and further expanded by H. Hurd [9] and W. Gardner [10].Since then, numerous scientific articles and books have addressed PC processes analysis.
In the literature, one can find various theoretical models that exhibit the PC property.In time series analysis, the most common is the periodic autoregressive moving average model (PARMA), which is a natural extension of the widely utilized autoregressive moving average (ARMA) model [11].In PARMA scenario, the constant coefficients in ARMA time series are replaced by periodic deterministic functions.In consequence, the periodicity of such models manifests itself in the ACVF.ARMA models find applications in various areas of interest, including engineering sciences, where they describe a real system based on measured data from mechanical systems [12][13][14].However, if the mechanical system operates under time-varying conditions, the coefficients of the ARMA model also fluctuate over time, displaying variability in some characteristics of the corresponding data.If such changes are periodic, it results in periodicity of data characteristics.Consequently, within mechanical systems, PARMA models (particularly PAR models) have been proposed as an optimal approach for systems experiencing cyclic variations in load or speed, such as the bucket wheel excavator used in the raw material industry [15].The application of PARMA and PC models in mechanical systems (particularly in condition monitoring) is the main motivation of our work.In this area, one of the key problems is the local damage detection in gears/bearings.Viewed through a signal processing perspective, damage detection can be characterized as the identification of impulsive signals (resulting from localized faults) and cyclic signals (stemming from the rotation of elements) within noisy observations.Thus, the identification of local damage can be considered as the detection of hidden periodicity.
In the classical analysis of PC processes, different approaches were used for the detection of hidden periodicity [16].However, the most common one is the analysis of the ACVF of the underlying signal in various domains (time, time-frequency, and bi-frequency domains) [17,18].The classical tool in such a case is the two-dimensional spectral density, a bi-frequency statistic, with its estimator called sample coherence.It is strongly related to the commonly used cyclic spectral coherence statistic, which serves as a tool for the identification of periodic characteristics of a signal [19].Methods for detection of hidden periodicity based on such maps are widely discussed in the literature.For example, the authors of [20,21] used new algorithms for cyclic spectral approaches in mechanical systems, while the study presented in [22] proposed a novel method for cyclic spectrum analysis improved by machine learning techniques.Several cyclostationary indicators (in the context of gear fault monitoring) that also make use of the ACVF were proposed in [23], see also [24,25].
Sample coherence, analyzed for example in [18,26], takes advantage of the property of PC models for which the support of the two-dimensional spectral density is a set of parallel, evenly spaced diagonal lines.It was also utilized for other statistics used to identify the hidden periodicity of PC time series, such as coherent and incoherent statistics [5,18].They were proposed in [27] as functions represented in the frequency domain to determine the presence of periodic correlation in given data.Both are one-dimensional representations of the two-dimensional sample coherence, which convey the information about potential periodicity in a way that is easier to interpret.When a periodicity is present in the analyzed signal, the correspondingly periodic peaks should be present in the plots of these statistics.
Although the methods mentioned above for the PC behavior identification appear to be effective for classical cases, in practical applications, it is frequently observed that signals display PC behavior (or, in general, hidden periodicity) but are also disturbed by additional noise.Consequently, detecting hidden periodicity becomes considerably more challenging compared to situations where a pure model (i.e., without the additive noise) is considered.Indeed, in industrial settings, the recorded signals are almost invariably corrupted by noise that can stem from disruptions in the measuring systems or other external factors that influence the observations.The problem of PC models with additive noise was discussed in the literature, where various techniques were proposed for the estimation of such models.For example, in [28,29] the authors examined the PAR models with additive outliers, and in [30,31], the noise-corrupted PAR system was discussed for the general distribution of noise with finite-variance distribution.The special case of noise-corrupted PARMA models (i.e., ARMA model) has been extensively discussed for different types of additive noise; see, e.g., [32][33][34].
In this paper, we discuss the problem of hidden periodicity detection when the model under consideration exhibits PC property but is disturbed by non-Gaussian additive noise.Similar to the classical approach (i.e., when the additive noise does not occur), the methodology is based on the sample coherence; however, in the algorithm we propose to use the robust version of the discrete Fourier transform based on the Huber function-based M-estimation.Using the robust version of the sample coherence, we propose robust counterparts of coherent and incoherent statistics and apply them to the hidden periodicity identification.Relatively simple modification of the classical approaches allows us to effectively apply them for noise-corrupted models.
Let us note that the problem of detecting hidden periodicity can be formulated using the language of hypothesis testing.In classical methodology based on Gaussian models, this approach seems to be natural.However, our situation differs from the classical approach because our goal is to demonstrate that detecting hidden periodicity is more complex when general additive non-Gaussian noise is present in the data.In hypothesis testing, we first need to define the test statistic.In our case, when considering the sample coherence-based statistic, it can serve as a basis for hypothesis testing.However, the next step involves identifying the test statistic's distribution under the null hypothesis.This is the main challenge because, under the null hypothesis, we do not assume any specific distribution of the data, nor do we assume any specific PC model.Therefore, formulating the problem in the language of hypothesis testing for such a general situation is difficult.In our previous paper (see [30]), we discussed the problem of hidden periodicity detection using hypothesis testing for the finite-variance case, where we assumed specific PC models.Here, the situation is more complicated, and the graphical methods proposed in this paper can help to identify hidden periodicity behavior.
In our analysis, we discuss two exemplary PC models.The first model (called later PC model 1) is a Gaussian noise with added periodic function, while the second one (called later PC model 2) is a Gaussian PAR model.Both models were explored in the literature as common systems for PC phenomena description in various areas.We also consider two types of additive noise.The first one (called AN model 1) is a sequence of additive outliers, while the second one (called AN model 2) is an autoregressive model with innovations coming from a general class of α−stable distribution [35].Let us note that the interdependence contained in AN model 2 makes the detection of hidden periodicity even more difficult.We recall that the α−stable distribution, in general, has infinite variance (except in the Gaussian case) and belongs to the so-called heavy-tailed class of distributions.This class has found many applications, including the condition monitoring area [36,37].Both considered PC models combined with two types of additive noise serve the stochastic systems with non-Gaussian distribution; however, in the case of AN model 1, this distribution has finite variance, while in the second case, it may have infinite-variance distribution.We highlight that, in real applications, additive noise may have a much wider sense than just a few outliers, and thus, we also consider the non-Gaussian α−stable additional noise in the considered models.The novelty of our approach lies in introducing robust algorithms for coherent and incoherent statistics through the application of the discrete Fourier transform based on Huber functionbased M-estimation and proposing a new technique for period estimation based on this methodology.Furthermore, we propose a robust version of the measure of fitness (MoF) statistic, which was originally discussed in [5].To our knowledge, incorporating a robust algorithm for discrete Fourier transform calculation has not been considered in the literature for coherent and incoherent statistics algorithms.Additionally, the technique for period estimation based on these robust algorithms is a novel proposition in the context of non-Gaussian additive noise.
We note that the concept of PC models with non-Gaussian distribution is not new; however, the robust coherent/incoherent statistics (obtained by using Huber functionbased M-estimation) is a new approach in this area.In the literature, there are known approaches for the analysis of signals exhibiting both hidden periodicity (cyclostationarity) and non-Gaussian behavior.The methods based on p-th-order cyclostationarity were proposed, e.g., in [38] for the direction of arrival estimation and in [39] for the time difference of arrival estimation.In the latter paper, the problem analyzed was also addressed using the fractional lower-order cyclostationarity; see also [40].This concept was also used in [41] for the joint estimation of time difference of arrival and frequency difference of arrival.In [42], phased fractional lower-order cyclic moments were considered.Another measure, called cyclic correntropy, was proposed in [43] and further analyzed, e.g., in [44,45].This approach was also used to propose robust methods for various signal processing problems, e.g., automatic modulation classification [46], cyclic array beamforming [47], and direction of arrival estimation [48].Methods based on the hyperbolic tangent cyclic correlation [49,50] and generalized cyclic correlation [51] were also considered.In [52], nonparametric cyclic detectors based on sign and rank statistics were proposed.
In this paper, by extensive Monte Carlo simulation studies for both PC models and both types of additive noise with different levels of non-Gaussianity, we demonstrate the efficiency of the proposed methodology based on robust MoF, coherent and incoherent statistics for hidden periodicity detection.Comparative analysis (classical versus robust versions of the analyzed statistics) clearly confirms the superiority of the proposed approaches in the considered cases.The simulation results are supported by analysis of real datasets from two areas.All considered time series have been previously discussed in the literature, also in the context of PC models with additive noise.The first dataset was examined in [53].It describes air pollution data, and in the mentioned bibliography position, the authors proposed to apply the PAR model with additive outliers (in our case it is PC model 2+AN model 1).In the current paper, based on the proposed robust statistics, we confirm the hidden periodicity property of the data and the superiority of the robust statistics over classical approaches.The other two datasets come from the condition monitoring area and describe the vibration signal from the crushing machine (operating in the raw material industry) and the acoustic signal from belt conveyor idlers, respectively.In our case, we apply the proposed robust statistics (coherent and incoherent) for time-frequency representations of the signals and, in consequence, we propose their new bi-frequency representation.We assume that both signals correspond to PC model 1+AN model 2, as they were analyzed in the literature [54,55].By the performed analysis, we confirm the hidden periodicity of the datasets and the superiority of the robust statistics over classical approaches in the considered context.The new bifrequency maps based on robust coherent and incoherent statistics can serve as efficient tools for local damage identification and may be effectively used for signals disturbed by non-Gaussian noise.
The rest of the paper is organized as follows.In Sect.2, we describe the considered PC models and the analyzed types of additive noise.In Sect.3, we recall the sample coherence and three statistics based on it-coherent, incoherent, and MoF statistics.We also introduce robust versions of the discussed statistics, which are based on the robust algorithm for the Fourier transform calculation.Here, we also describe a new algorithm for period estimation, which is based on robust statistics proposed in this paper.In Sect.4, we present the Monte Carlo simulation study, where for different levels of non-Gaussianity of additive noise we demonstrate the superiority of the proposed approaches over classical methods.In Sects.5 and 6, we present the analysis of real datasets.The last section concludes the paper.

Periodically correlated time series with additive noise
In this paper, we consider the following model, called periodically correlated model with additive noise where {X t } is a periodically correlated (PC) time series (random sequence) with period T ∈ N * and {Z t } is the additive noise (AN).We assume that {Z t } is a stationary time series with a non-Gaussian distribution that is independent of {X t } .First, we recall the definition of PC time series.

Definition 1
The time series {X t } is periodically correlated (or second-order cyclosta- tionary) if its mean and autocovariance functions are periodic in t with the same period T ∈ N * [18] (1) where E(•) is the expected value, cov(•, •) is the covariance, and there are no smaller val- ues of T ∈ N * for which both conditions in Eq. ( 2) hold.
As a particular case of the periodically correlated time series, we consider the model for which only the mean or the autocovariance function is periodic while the other characteristic (i.e., autocovariance function and mean, respectively) is independent of time t (i.e., the mean function is constant in t, and the autocovariance function depends only on lag h).In such a case, the time series still fulfills Definition 1.Note that a function independent of time t is periodic with any period T ∈ N * .

Examples of PC time series
In this paper, we consider two examples of time series that exhibit the PC property given in Definition 1.
PC model 1.As the first model of PC time series {X t } , we consider the sum of a weakly stationary time series and the T-periodic deterministic function f : Z → R (i.e., f (t + T ) = f (t) ).Then, the random sequence {X t } defined as is a PC time series, see [56].The most typical example of a weakly stationary random sequence {ξ t } used in Eq. ( 3) is the sequence of Gaussian N (0, 1) independent identically distributed (i.i.d.) random variables.This is also the case considered in this paper.
PC model 2. The second considered example of a PC model is the periodic autoregressive (PAR) time series considered as the periodic extension of the well-known autoregressive (AR) model that assumes constant coefficients.The PAR(p) time series is defined as follows [18] where the {ξ t } is a sequence of i.i.d.random variables.In this paper, we assume that for each t, ξ t ∼ N (0, 1) .The parameter sequences {φ i (t), i = 1, . . . ,p} are periodic with the same period T ∈ N * with respect to t.

Examples of additive noise
In this paper, we consider two types of additive noise.In both considered cases, we assume that the time series {Z t } in Eq. ( 1) is stationary and has a non-Gaussian distribu- tion.However, the AN models considered exhibit different properties.
AN model 1.The first considered model of AN assumes that {Z t } in Eq. ( 1) constitutes a sequence of i.i.d.random variables such that for each t we have where {K t } is an i.i.d.sequence of the following distribution with P denoting the probability.We assume A t and K t are independent for each t.Moreover, {A t } also constitutes an i.i.d.sequence and for each t, A t is a random variable from the uniform distribution additive outliers sequence.For each t, Z t has finite variance for any set of parameters, even though it exhibits impulsive behavior.

AN model 2.
As the second AN model, we consider the AR(p) time series with α− stable distribution.Precisely, in this case, we assume that {Z t } in Eq. ( 1) satisfies the following equation [57,58] where the {ε t } is a sequence of i.i.d.random variables with symmetric α−stable distribu- tion defined via the characteristic function [59,60] In the above, 0 < α ≤ 2 represents the so-called stability index and σ > 0 is the scale parameter.Let us recall that the symmetric α-stable distribution is considered as the classical example of the heavy-tailed (power-law) distributions.Moreover, it is an extension of the Gaussian distribution, as for α = 2 it reduces to this classical distribution.The stability index is responsible for the heaviness of this distribution's tail.The smaller α , the higher the probability of large values.For α < 2 , the variance of the α-stable distri- bution is infinite.
Let us note that, in contrast with the AN model 1, in this case the sequence {Z t } does not constitute a sample of independent random variables as long as at least one parameter of the AR model is nonzero.In general, {Z t } exhibits the so-called short- range dependence [58].The distribution of Z t for each t belongs to the non-Gaussian family of distributions (here α−stable).In the case of α = 2 , the α−stable AR model defined in Eq. ( 7) reduces to the classical AR model with Gaussian innovations; see [11].Finally, when p = 0 , the discussed AN model is just a sequence of i.i.d.random variables with strictly α−stable (when α < 2 ) or Gaussian (when α = 2 ) distribution.
The following setups of presented PC models (with period T = 8 ) and additive noise models are considered The sample trajectories (of length N = 720 ) of the PC1 and PC2 models and their combinations with both the AN1 and AN2 sequences are presented in Fig. 1.One can see that the periodic behavior cannot be detected by visual inspection even in PC1 and PC2 samples, that is, without any additive noise.However, in these cases, it would be easily identified using classical tools, such as those presented in Sect.3.1.On the other hand, when additive noise is included, it significantly corrupts the behavior originally observed.In trajectories with AN1, only several observations are changed in comparison with the original PC time series, but they are significantly outlying, which strongly influences results of standard methods negatively.In the AN2 case, all (7) observations are changed in comparison with the original PC time series, and some peaks can be observed as well-hence, it can be considered a more challenging case.

Classical methods for periodicity detection
In this paper, we focus on periodicity detection algorithms in the bi-frequency domain.
All considered methods are based on the sample coherence statistic, which is a smoothed estimator of the normalized two-dimensional spectral density function [5].This approach is one of the most commonly used for the detection of hidden periodicity, especially for the models where the autocovariance function is periodic with period higher than 1.However, our goal is to highlight the universality of the introduced methodology; therefore, we also apply the sample coherence statistic for the models where only the mean function is periodic (while the autocovariance function is independent of time t).For a time series X = [X 0 , . . ., X N −1 ] , the sample coherence statistic is defined as where B is the smoothness coefficient, and I(ω j ) = N −1 t=0 X t exp(−2π iω j t) is the dis- crete Fourier transform for standard frequencies ω j = j/N (measured in cycles/sample), j = 0, . . ., N − 1 .The usefulness of this statistic comes from the fact that for PC models, the support of the two-dimensional spectral density function is a set of parallel, evenly spaced diagonal lines.The two-dimensional plot of sample coherence may serve as a visual inspection of periodicity presence.
The sample coherence is most often only a base for other tools, which produce results that are easier to interpret and apply in automated procedures.In [27], the coherent and incoherent statistics were proposed.Both can be considered as one-dimensional representations of the originally two-dimensional sample coherence.For 0 < d < N , the coherent statistic is defined as (9)  The choice of this parameter is discussed in [18,27].Note that δ I (d) is a generalization of δ C (d) and simplifies to the latter for B = N .Another method of periodicity detection is the measure of fitness (MoF) statistic [5].It is a bootstrapbased approach that, in some cases, might be more efficient than coherent and incoherent statistics.The MoF is defined as for a specified value of B, with where c is the estimated critical value (for a given significance level) obtained using the moving block bootstrap (MBB) method [61].All three presented statistics produce values from the interval [0,1].In the above formulation, they are functions of the d = |l − k| parameter, that is, in the domain of frequency indices.Then, the peaks at d = d T , 2d T , 3d T , . . .values indicate a perio- dicity with period T = N /d T .However, these statistics are often considered in the frequency domain, i.e., as functions of ω d = d/N .This argument can also be further adjusted, for example, to include the sampling frequency of the analyzed time series.

Robust methods for periodicity detection
Although the statistics presented above serve as efficient tools for periodicity detection in non-impulsive time series (in particular Gaussian), all of them (as will also be demonstrated later) may be inefficient when applied to data with non-Gaussian, impulsive behavior.The reason for this is the sensitivity of the discrete Fourier transform I(ω j ) to outliers [62].To solve this problem, we consider here a robust approach for I(ω j ) calculation proposed in [63].It is based on the fact that the Fourier trans- form calculation can be formulated as a linear regression problem in the following way where the coefficients βj = [ β1,j , β2,j ] ′ are obtained by minimizing the least square cost function (10) for C r = [cos(2π rω j ), sin(2πrω j )] ′ .The terms X r − C ′ r β j , r = 0, . . ., N − 1 , are later called residuals.The regression in Eq. ( 15) is applied only for j = 1, . . ., ⌊(N − 1)/2⌋ .The results for larger frequencies can be obtained using the property I(ω j ) = I(1 − ω j ) , as we only consider real time series X .Moreover, for ω j = 0 and ω j = 1/2 , when sin(2πrω j ) = 0 for each r ∈ N , we set β2,j = 0 and apply the regression only to find the β1,j term.For these ω j , we then set I(ω j ) = N β1,j .
The least squares method is very inefficient when outliers are present in the data; however, there are various approaches for robust linear regression which aim to overcome this drawback.In this paper, we use the idea of M-estimation [64], where the square cost function in Eq. ( 15) is replaced by a function ρ(•) that is less sensitive to outliers.In other words, we consider the following optimization problem which is equivalent to solving where ψ(•) = ρ ′ (•) , and s = MAD/0.6745 is the selected robust scale estimate with median absolute deviation (MAD) of the residuals' vector in the numerator.Let us note that ρ(•) and ψ(•) functions are applied to the residuals standardized using s.
There are many ρ(•) functions with robust properties, e.g., Huber, Tukey's biweight, and Andrews functions.In this paper, we use as ρ(•) the Huber function, which is defined as where θ > 0 is a tuning constant.One can see that this function is built upon square ( L 2 ) and absolute value ( L 1 ) loss functions, combining advantages of both, i.e., efficiency in standard cases (without outliers) and robustness, respectively.The trade-off between them is determined by the value of θ .The lower it is, the more robust (but less efficient in standard cases) the method is.Usually, we assume θ = 1.345 , for which the Huber regression is 95% as efficient as the least squares method in the Gaussian case.In practice, M-regression is usually performed using the iteratively reweighted least squares (IRLS) method [65].In this paper, we used the robustfit function in MATLAB.The robust Fourier transform Ĩ(ω j ) can then be derived using Eq. ( 14), i.e., we set where βj are now calculated from Eq. ( 16).Again, we apply this regression only for j = 1, . . ., ⌊(N − 1)/2⌋ and use Ĩ(ω j ) = Ĩ(1 − ω j ) to obtain the results for larger (15) frequencies.For ω j = 0 and ω j = 1/2 , we set β2,j = 0 and look only for the β1,j term, finally setting Ĩ(ω j ) = N β1,j .
Using Ĩ(ω j ) , we are able to modify the periodicity detection algorithms described in Sect.3.1 to make them more robust and suitable for non-Gaussian time series.The only difference in comparison with their classical versions is that all are based on the robust sample coherence defined as It can be seen that the above formula is a modification of Eq. ( 9) with all terms I(ω j ) replaced by the corresponding Ĩ(ω j ) values.
Based on robust sample coherence, we define robust coherent and incoherent statistics, respectively, as Analogously, we define the robust MoF as where c is the estimated critical value (for a given significance level) obtained using the MBB method.

Method for period estimation
As mentioned above, for all statistics presented (in both standard and robust versions) the peaks at d T , 2d T , 3d T , . . .indicate periodic behavior with period T = N /d T .For all other d, we expect the obtained values to be significantly lower.This property of the considered methods allows us to propose a simple criterion for period selection in an automatic manner.Let us assume that for a time series X one of the above statistics (in standard or robust version) δ(d) was calculated for d = 1, . . ., d max .Note that for any potential period T * (for simplicity, let us assume that T * divides N), the expected placements of the aforementioned peaks are known.Then, we define the period selection criterion as Then, out of all considered potential periods T * , we select the one for which PSC(T * ) takes the largest value.If there is a tie, we choose the smallest potential period.Let us note that the methodology for period estimation presented above is universal and can utilize both classical and robust statistics discussed in Sects.3.1 and 3.2. (20 4 Monte Carlo simulations

Periodicity detection
This part of the simulation study is devoted to the efficiency analysis of the periodicity detection methods described in Sects.3.1 and 3.2.Throughout entire Sect.4, for both PC models and both AN models considered, we assume the same settings as for the exemplary trajectories presented in Fig. 1.Let us recall them below • PC model 1 (later denoted as PC1): Therefore, as follows from the construction of PC1 and PC2, in all cases considered we have the period T = 8.

PC model 1
Let us first consider the PC1+AN1 and PC1+AN2 cases.For the trajectories of these models presented in Fig. 1 (recall their length N = 720 and period T = 8 , thus d T = 90 ), we calculate each of the six considered statistics for d = 1, . . ., d max = 360 .In the incoherent and MoF methods (in both versions), for all experiments presented in this paper, we set B = 60 .Moreover, for MoF methods, we use 100 bootstrap replicates and the 1% signifi- cance level.The results are illustrated in Fig. 2. In these plots, red dashed lines are included for d ∈ D = {90, 180, 270, 360} (called later cyclic d).These are the expected placements of statistics' peaks (relative to values obtained for other d) up to d max , for the periodicity of T = 8 to be correctly identified.Moreover, to quantify its efficiency, for each statistic (denoted in the formula below as δ(•) ), we calculate the following indicator (24) τ = d∈D δ(d) Most of all, one can see a clear difference between the results of the standard and robust version of each statistic.All classical statistics behave in a way that does not indicate the actual periodicity of the underlying time series.On the other hand, in their robust versions, a clearly visible peak for d = 90 is present, which is enough to correctly identify the periodicity with T = 8 .The highest efficiency indicator τ was obtained for the robust coherent statistic, for which the aforementioned peak is the most significant (relatively to other values observed in its plot).The results for the PC1+AN2 model confirm that this case is more challenging than the previous one.For all robust statistics, the periodicity is less distinguishable than before, as the peaks in d = 90 are now relatively less significant.However, robust coherent and robust MoF statistics, for which the highest values of τ were obtained, still performed acceptably.On the other hand, all standard statistics again did not detect periodicity.In fact, let us note that for both PC1+AN1 and PC1+AN2 cases, the obtained values of τ were always higher for robust methods than for the corresponding standard approaches.
Next, let us analyze the efficiency of the proposed methodology for different parameters of both considered additive noise distributions, using the proposed performance indicator τ and Monte Carlo simulations.In the following comparisons, we exclude both standard and robust MoF methods because of their computational complexity.Figure 3 presents the values of efficiency indicator τ averaged over 1000 simulated trajectories from PC1 model with given additive noise.In the left column, it is the AN model 1 with q = 0.005 for different values of a (top panel) and with a = 60 for differ- ent values of q (bottom panel).Analogously, in the right column, it is the AN model 2 with φ 1 = 0.2 (AR(1)), with σ = 1 for different values of α (top panel) and with α = 1.8 for different values of σ (bottom panel).First, let us analyze the results obtained for AN model 1.Most of all, let us note that for more difficult cases, related to a larger magnitude and frequency of outliers in a time series (i.e., when a and q increases), both considered robust methods maintain their efficiency, with a large advantage of robust coherent statistic.On the other hand, for standard methods, the efficiency tends to decline.These results clearly confirm what could also be seen before, the sensitivity to outliers of non-robust methods, and the lack of this drawback for robust approaches.For AN model 2, the efficiency of both robust methods also tends to decrease more slowly for more extreme cases (that is, for lower α and larger σ ).The observed results again show the significant advantage of the robust coherent method over all other approaches.

PC model 2
Let us now perform the same analysis, but for PC2 case as the underlying periodically correlated model.Figure 4 presents the results for all six considered statistics applied to PC2+AN1 and PC2+AN2 samples illustrated in Fig. 1.Moreover, in Tab. 2, we show the values of the corresponding efficiency indicator τ .As before, for the case with AN1, the advantage of the proposed robust methods is very significant.For all of them, the observed peaks clearly indicated periodicity with T = 8 , whereas their standard versions failed, showing once again their strong sensitivity to outliers in a time series.Non-robust methods did not perform well neither for the model with AN2.Moreover, let us note the inefficiency of the robust coherent method.However, both robust incoherent and robust MoF approaches again succeed in the periodicity detection, with their peaks placed in an expected way.For both PC2+AN1 and PC2+AN2 cases, the best efficiency indicator τ value was obtained for the robust MoF method.Again, for both cases, each robust method yielded a larger value of τ than its standard counterpart.
Let us now further investigate the performance of the analyzed methods using Monte Carlo simulations.The plots of the averaged efficiency indicator for the PC2 model and different parameters of additive noise (analogous to Fig. 3 for PC1) are presented in Fig. 5.For the AN model 1, again, one can see that both robust methods are insensitive to outliers, and hence, they are much more efficient than standard techniques.This is also true for AN model 2, where for all parameters considered (except for the Gaussian case α = 2 ) robust approaches outperform the classical ones.Let us note that here, as opposed to the PC1 model, the robust incoherent statistic was more efficient than the robust coherent statistic.

Period estimation
In this part, we analyze the efficiency of the period selection criterion presented in Sect.3.3 based on different periodicity detection statistics δ(d) .Here, we evaluate only coherent and incoherent statistics with their robust versions (again omitting both MoFbased methods due to their computational complexity).

PC model 1
Let us first consider the PC1+AN1 case.We simulate 1000 trajectories of length N = 720 of this model.Then, for each sample, we calculate each statistic δ(d) for d = 1, . . ., d max = 360 and then the PSC(T * ) criterion based on it for the period can- didates T * = 2, 3, 4, 5, 6, 8, 10, 12, 15, 16 (recall that the true period is T = 8 ).Boxplots of the criterion values obtained for all generated trajectories, for each δ(d) statistic and period candidate T * , are presented in the left panel of Fig. 6.The results of the analogous procedure for PC1+AN2 case are illustrated in the right panel.Here, for each method, we expect the criterion values for T * = T = 8 to be significantly higher than for the other candidates.One can see that such behavior is clearly visible for both robust methods, and the proposed criterion based on them can indeed identify the correct period.
Although the results for T * = 16 are at the same level, let us note that it is double the correct period and, due to the construction of PSC(T * ) , it will always have at least the same value of this criterion (hence, in case of ties, we always select the smallest period candidate).For both the classical coherent and the incoherent statistic, the criterion is less efficient, since the results for T * = 8 are not significantly larger than the others (especially in the PC1+AN1 case).Similar to the analysis of periodicity detection statistics, let us now analyze the performance of the proposed period selection criterion for different parameters of both AN distributions.For that purpose, we conduct the following experiment.For the PC1 model with given additive noise, we simulate 1000 trajectories of length N = 720 .For each sample, we calculate a given statistic δ(d) for d = 1, . . ., d max = 360 and the PSC(T * ) based on it for T * = 2, 3, 4, 5, 6, 8, 10, 12, 15, 16 to identify a period for this tra- jectory.In the end, we calculate the percentage of cases (out of 1000 simulated trajectories), for which the estimated period is equal to the true one T = 8 .As δ(d) , we again consider only coherent and incoherent statistics (in classical and robust versions).The results for different parameters of AN model 1 and AN model 2 (arranged in the same way as, for example, in Fig. 3) are presented in Fig. 7. Most of all, one can see that for all cases considered with AN model 1 (with different a and q), the robust statistics-based criteria are almost perfectly efficient.A similar behavior can also be observed for AN model 2 with different values of α , with a small decrease for higher values of σ .This confirms that the proposed period selection criterion combined with robust periodicity detection methods is able to serve its purpose.On the other hand, the analogous results for classical statistics again prove their inefficiency for considered models-one can see a strong decline of the success percentages as the additive noise becomes more significant.

PC model 2
Let us now conduct the same analysis for the PC2 model.The boxplots of the PSC(T * ) criterion values calculated for the PC2+AN1 and PC2+AN2 models are presented in Fig. 8.One can see that these cases (in particular the latter) turned out to be more challenging than the PC1-based ones.For both classical statistics, the results for T * = T = 8 again do not stand out from the others.Moreover, for the robust coherent statistic, the results for T * = 4, 12 are comparable to those for T * = 8, 16 which could lead to a wrong period estimation.The desired behavior can only be observed for the robust incoherent statistic, in particular for the PC2+AN1 model.In the case of the presence of AN2, the results for the correct period are still visibly greater than for others, although the difference in comparison to other even period candidates is not large.The percentages of correctly identified cases for the PC2 model with different additive noise parameters are presented in Fig. 9.For AN model 1, again one can see that the efficiency of both robust methods does not change significantly for different parameters a and q, although the robust incoherent statistic is always much better with almost perfect effectiveness.It is also the best method for all the setups considered of AN model 2 (with a slight exception for the Gaussian case α = 2 ), even though its performance decreases for more extreme cases (especially for larger σ ).As before, both classical methods are getting significantly worse for larger levels of additive noise, and in general, the advantage of robust approaches is very clear.

Application for air pollution data
Let us now present the application of the proposed periodicity detection methodology to real data from the air quality area.The time series analyzed are the daily mean pollutant particulate matter with diameter smaller than 10 mm (PM 10 ), measured hourly in µg/m 3 and collected at the station located in the Great Vitória Region GVR-ES, Brazil, at the Automatic Air Quality Monitoring Network (RAMQAr).A complete description of the data can be found, for example, in [53].In this article, it was claimed that the described phenomenon can be modeled by the PAR time series with additive outliers, which in our case corresponds to the PC2+AN1 combination considered in this paper.This hypothesis was confirmed in [30] using the proposed methodology for testing the presence of additive noise.In both mentioned positions, the period T = 7 was assumed a priori because of a weekly rhythm of the analyzed phenomenon.In this paper, we apply the proposed methodology to validate this assumption and compare the results of standard/robust approaches.We consider two daily mean PM 10 time series: from Carapina and Vitória (center) monitoring stations, both gathered from January 1, 2018, to May 5, 2019 and consisting of 490 observations.They are presented in Fig. 10.In both series, one can identify some outlying observations; however, there are more of them in the Carapina dataset.
For both considered datasets, we apply the proposed periodicity detection statistics in a similar manner as for the exemplary trajectories in Sect.4.1 (see Figs. 2 and 4), that is, all six considered statistics are applied to each time series, taking d = 1, . . ., d max = 245 .The results are presented in Fig. 11.The red dashed lines now correspond to the expected placements of the statistics' peaks for the period T = 7 , which are d ∈ D = {70, 140, 210} .
In the case of the Carapina dataset, one can see that both versions of coherent statistics returned a peak for d = 70 as expected; however, in the robust case it is a little bit more significant in comparison with other obtained values of a statistic (in particular, to the peak obtained around d = 100 ).On the other hand, in both incoherent statistics, there are no peaks that would indicate a periodicity.Both MoF statistics returned a peak in d = 70 , but some other peaks are also present.
The results for the Vitória (center) dataset show that, in this case, it was easier to identify the periodicity.Again, coherent and robust coherent statistics returned a peak in d = 70 .For both incoherent statistics, large values were obtained for d = 70, 140 , Fig. 10 Two analyzed air pollution datasets-daily mean PM 10 measured in Carapina and Vitória (center) stations although they are relatively less significant than the peaks of both coherent statistics.In the case of both MoF methods, for all d = 70, 140, 210 the obtained values are large, but again, other peaks are present as well.
For both considered cases, the efficiency indicator τ (see Eq. ( 24)) was calculated.The results are presented in Tab. 3. One can see that, in general, the indicator values for robust methods are larger than for their non-robust counterparts.The exceptions are incoherent statistics for the Carapina dataset, where both standard and robust versions failed, and MoF-based statistics for Vitória (center) data.In the latter case, a small advantage of the standard method may be caused by the fact that, as presented, for this dataset the periodicity detection was generally easier, and the influence of outlying values was not significant.

Application for local damage detection problem
In this section, we present the application of the proposed methodology to real datasets from the condition monitoring area.We analyze two datasets corresponding to vibration signals from a crushing machine and acoustic signal from belt conveyor idlers, which may be considered as bases for the identification of local damage in rolling element bearing.There are many approaches for fault modeling in bearings [66,67].For simplicity, it is assumed that a signal containing information about damage is a sum of deterministic periodic function (local impulses) and random noise.Let us note that this case corresponds to PC model 1+AN model 2 considered in this paper.
For both analyzed signals, we will analyze the performance of the proposed methodology when periodicity is present in the signal (case of local fault), as well as when it is not (healthy case).In the presented analysis, we again consider only coherent and incoherent statistics ( B = 60 ) in classical and robust versions, excluding both MoF methods.Although the considered methods can be used in a straightforward way for analyzed signals in the time domain, here we present an alternative approach where they are applied to signal represented in the time-frequency domain (spectrogram).This idea is often used in local damage detection methods because of the very complex structure of the analyzed signals.The local impulses corresponding produced by a damage are often identifiable only in a specific frequency band (called the informative frequency band, IFB).Hence, here we also use the time-frequency representation.Such a procedure is also much more time efficient, because the statistics are not applied to a single, very long data vector, but to a number of much shorter samples, which significantly reduces the computational complexity.Moreover, as a result of this approach, we obtain a bi-frequency map, which allows us to easily identify the periodicity as well as the informative frequency band (in a similar manner as, e.g., spectral coherence maps [19]).
For each considered statistic δ(d) , the procedure for creating a spectrogram-based bi- frequency map is as follows.First, for the analyzed signal X = [X 0 , . . ., X N * −1 ] we calcu- late the spectrogram S(t, f) where STFT(t, f ) is the short-time Fourier transform given by the formula with w(t − n) being a shifted window of length N w .Then, for each frequency f, we take the corresponding spectrogram row S f = [S(t 1 , f ), . . ., S(t N , f )] (its length N depends on selected spectrogram parameters), and for selected ǫ min and ǫ max , we calculate for each ǫ min ≤ ǫ ≤ ǫ max (modulation frequency) the statistic δ(d) setting d = ǫN * /f s , where f s is the sampling frequency.As a result, we obtain a bi-frequency map denoted as �(f , ǫ) .To compare all obtained maps with the same scale [0,1], we rescale each map by dividing each of its values by the maximum value of this map.

Vibration signal from a crushing machine
The analyzed vibration signal from a crushing machine is a two-second sample from the dataset considered in [54].The system used for measurement consists of Endevco accelerometers (vibration), BruelKjaer Laser probe (shaft speed was measured), NI DAQ card, Labview SignalExpress and laptop.The vibration signals were acquired in horizontal and vertical directions.The analyzed signal and its spectrogram (constructed using nfft = 512 frequency points and Hann window of length N w = 128 and overlap 110) are illustrated in Fig. 12.This case is further referred to as Dataset 1a.In the plot of the signal, one can see a significantly outlying impulse around 1.5s, which might cause a challenge for standard, non-robust methods.The length of this signal is N * = 50000 , and the sampling frequency is f s = 25000 Hz.It is important to note that this signal was collected from a healthy crushing machine, which means that no fault periodicity is expected in the results.
Using the presented spectrogram, let us construct the analyzed bi-frequency maps for each of four considered statistics, taking ǫ min = 3 Hz and ǫ max = 100 Hz.The results are (25) presented in Fig. 13.Most of all, one can see a clear difference between standard coherent/incoherent statistics and their robust versions.In the former, unexpectedly, very significant values were obtained for f > 5000 Hz, whereas in the latter, the constructed maps are much more homogeneous, which is a more desired result in this case.
As mentioned, Dataset 1a was gathered from a healthy machine; hence, no periodicity is present there.To imitate a damage of the machine, let us artificially introduce the periodicity to the original signal, by adding a series of cyclic impulses of a form of decaying harmonic oscillation with amplitude A = 45 , informative frequency band f c = 3500 − 6500 Hz, and cyclic frequency f f = 30 Hz.The resulting signal (called Dataset 1b) and its spectrogram (constructed using the same parameters as in Dataset 1a) are presented in Fig. 14.One can see that the introduced cyclic impulses are covered by the background, non-cyclic noise-in particular, by the mentioned outlying impulse around 1.5s.
The considered bi-frequency maps calculated for Dataset 1b are presented in Fig. 15.As now the periodicity is present, we expect all maps to detect it by returning significantly high values for ǫ = 30, 60, 90 (true cyclic frequency and its multiples), in the 3500 < f < 6500 band, and low values everywhere else.The results obtained show that this is the case only for both considered robust methods.For classical coherent/incoherent statistics, the significantly large disturbances for f > 5000 (present also in the results for Dataset 1a) cause that the periodicity cannot be detected.Therefore, the advantage of the proposed robust techniques is clear.In particular, the periodicity is easily identifiable in the robust coherent statistic-based map.
To quantify the performance in periodicity detection, we calculate the following indicator for each constructed map where the larger the value, the more efficient the method is.This indicator reflects the desired behavior of each map in a similar way as the previously considered τ (Eq.( 24)).Let us note that again it is the ratio of the sum of values obtained for arguments corresponding to the correct periodicity detection (the selected f and ǫ range is explained (27)  Let us now analyze the signal from a damaged machine, where periodicity is naturally present.This signal (later referred to as Dataset 2b) and its spectrogram (constructed using the same parameters as for Dataset 2a) are illustrated in Fig. 18.Again, the signal length is N * = 96000 with the sampling frequency f s = 48000 Hz.Here, one can iden- tify some cyclic impulses; however, they are covered by two strong non-cyclic impulses around 0.6s and 1.8s.Their negative influence on the results of standard methods is significant, as shown in Fig. 19, where the maps for Dataset 2b are presented.For both nonrobust methods-based maps, the values obtained in some frequency bands are too large for the periodicity to be correctly detected.On the other hand, both robust methods much better recovered the information about periodicity.Although the obtained cyclic impulses are not as clearly observable as in Dataset 1b case, their visibility is sufficient to identify a periodicity with f c = 5.5 Hz.As for Dataset 1b, let us also consider the perfor- mance indicator for each map, now given by It is constructed analogously as τ 1 (�) (Eq.( 27)) considered for Dataset 1b.The only dif- ference is related to constraints in the numerator-we expect all maps to return large values for ǫ = 5.5 and its multiples, without specified range for f (since here the inform- ative frequency band is not known).As presented in Tab. 5, the proposed indicator confirms the advantage of robust methods which obtained better results than their nonrobust counterparts.

Conclusions
In this paper, we discussed the problem of hidden periodicity detection when the model under consideration exhibits PC property and is additionally disturbed by non-Gaussian distributed noise.We proposed a modification of the sample coherence algorithm by incorporating the robust version of the discrete Fourier transform.This procedure enabled the introduction of robust versions of coherent and incoherent statistics, as graphical methods for detecting hidden periodicity when the PC model is affected by additive non-Gaussian noise.Furthermore, we analyzed the robust MoF statistic, which also utilizes the robust sample coherence.Building on this approach, we introduced an algorithm for period estimation specifically tailored for the considered class of models.The main goal of this paper was to propose a general (28)   In our simulation study, we analyzed two types of PC models and two types of additive noise.They were chosen to encompass the properties observed in the real data examined.We demonstrated the effectiveness of our proposed approach across various levels of non-Gaussianity of the additive noise.The results are compared with those obtained from a classical approach that employs the standard sample coherence.The presented simulation studies clearly validated the rationale for employing robust techniques with the models considered.The most challenging scenario emerged with the Gaussian PAR model (i.e., PC model 2) combined with additive noise modeled as an α-stable autoregressive time series (AN model 2).Incorporating the heavy-tailed α -stable distribution, which may exhibit infinite variance, along with the interdependence of additive noise (i.e., autoregression), posed a significant challenge in detecting hidden periodicity.However, we have demonstrated the utility of our methodology even in this complex case.
As an illustration of the proposed approach, we analyzed three datasets from the fields of environmental engineering and condition monitoring.The first dataset pertained to air pollution, which has been previously examined in the literature, where the PAR model with additive outliers was proposed for the analysis of such data.Leveraging this assumption, we applied robust techniques for hidden periodicity detection and demonstrated their effectiveness in this context.The other two datasets involved vibration signals from a crushing machine and acoustic signals from belt conveyor idlers.Both signals have previously been studied in the context of signal-based local damage detection and were modeled using a PC model with non-Gaussian α-stable noise.From a sig- nal processing perspective, the problem of local damage detection can be approached by the detection of hidden periodicity.We analyzed these signals using time-frequency representations (spectrograms) and applied robust coherent and incoherent statistics to these representations.As a result, we proposed new robust bi-frequency maps that provide clear information about local faults, even in the presence of impulsive disturbances.
Although the results were presented using only specific real-world datasets, the proposed methodology can be applied to other areas of interest, as the issue of hidden periodicity and external factors inducing non-Gaussian behavior in data is prevalent across various domains.
In our future work, we intend to further explore this research within the field of condition monitoring.Specifically, we aim to leverage the newly developed bi-frequency maps to design statistics that can aid in predicting local faults or identifying them at an early stage.In the literature, such statistics are commonly referred to as health indices and play a crucial role in prognostics for estimating the remaining useful life of machinery.

Fig. 1
Fig. 1 Sample trajectories of PC1 (left column) and PC2 (right column) time series: without additive noise (top row), and with additive noise from AN1 (middle row) and AN2 (bottom row) models

Fig. 3 Fig. 4
Fig. 3 Average efficiency indicator values for PC1 and different parameters of AN models (left column: AN model 1, right column: AN model 2)

Fig. 6
Fig. 6 Boxplots of period selection criterion PSC(T * ) values based on the analyzed periodicity detection statistics for different T * , for PC1+AN1 and PC1+AN2 cases

Fig. 7
Fig. 7 Percentages of cases with correctly identified periods for PC1 and different parameters of AN models (left column: AN model 1, right column: AN model 2)

Fig. 8 Fig. 9
Fig. 8 Boxplots of period selection criterion PSC(T * ) values based on the analyzed periodicity detection statistics for different T * , for PC2+AN1 and PC2+AN2 cases

Fig. 12
Fig. 12 Analyzed real vibration signal from a crushing machine (Dataset 1a) and its spectrogram

Fig. 14
Fig. 14 Modified vibration signal with added periodic impulses (Dataset 1b) and its spectrogram

Fig. 17 Fig. 18
Fig. 17 Bi-frequency maps �(f , ǫ) based on the spectrogram and coherent/incoherent statistics (in standard and robust versions) for the Dataset 2a

Fig. 19 Table 5
Fig. 19 Bi-frequency maps �(f , ǫ) based on the spectrogram and coherent/incoherent statistics (in standard and robust versions) for the Dataset 2b

Table 1
Values of efficiency indicator τ for all analyzed periodicity detection statistics for PC1+AN1 and PC1+AN2 cases (with the best results marked in bold)

Table 2
Values of efficiency indicator τ for all analyzed periodicity detection statistics for PC2+AN1 and PC2+AN2 cases (with the best results marked in bold)

Table 3
Values of efficiency indicator τ for all analyzed periodicity detection statistics for both considered daily mean PM 10 datasets assuming T = 7 (with the best results marked in bold)