HOS network-based classification of power quality events via regression algorithms

Abstract

This work compares seven regression algorithms implemented in artificial neural networks (ANNs) supported by 14 power-quality features, which are based in higher-order statistics. Combining time and frequency domain estimators to deal with non-stationary measurement sequences, the final goal of the system is the implementation in the future smart grid to guarantee compatibility between all equipment connected. The principal results are based in spectral kurtosis measurements, which easily adapt to the impulsive nature of the power quality events. These results verify that the proposed technique is capable of offering interesting results for power quality (PQ) disturbance classification. The best results are obtained using radial basis networks, generalized regression, and multilayer perceptron, mainly due to the non-linear nature of data.

1 Introduction

With the consequent unstoppable increase of electronic equipments demanding electricity, consumers expect uninterrupted availability and quasi perfect power quality (PQ). For that reason, PQ is being the object of continuous interest by researchers and developers, due to its influence over and from the loads, and the recent potential inclusion in the modern smart grid (SG).

In this frame, an adequate PQ assures the necessary compatibility between all equipment connected to the grid [1]. The terms of this compatibility gather several aspects: sustainable power with low losses and high quality and security of supply and safety, being at the same time economically efficient, reliable, and resilience [2,3].

Certainly, the future SG would introduce transformative technologies to meet these design requirements, integrating intelligence into end-use devices as the key to satisfying the demand response. In parallel, the industry would design ways to incorporate automatic end-use-load participation into the model so that customers are not bothered by these programs and decisions, and so their lifestyles are not inconvenienced [4].

Due to the above arguments, the role of smart meters and sensors is the first being revised in the present and future SG. These automated meters (AM) used a two-way communicating infrastructure and centralized management, as well as new features such as the following: outage management, demand response, automatic load shedding, distribution automation, and the ability to enable and commute alternative energy sources.

Provided with this scenario, this research integrates artificial neural networks (ANNs) and advanced signal processing techniques based in higher-order statistics (HOS) in order to be implemented into an automated smart meter for PQ event detection and classification, within the frame of a SG with a high distribution penetration of renewable sources. Seven regression algorithms are tested based on an hybrid time-frequency battery of characteristics, specially designed to deal with non-stationary measurement time series.

The feature extraction stage from PQ disturbances is based on HOS, idea which has been proven to be efficient in several works. Indeed, since PQ events are sudden changes in the power line, the HOS are potentially useful to characterize each type of electrical anomaly both in time and frequency domains. As a novelty, the present research exploits the combination of time and frequency domain features, to deal with the inherent non-stationary associated to the electrical anomalies, with the goal to improve the performance of the ANN and to make feasible the integration in an smart meter.

Regarding backgrounds of HOS applications in this field, in the time domain, several notable works are worthy, e.g., Bollen et al. introduced new advanced statistical features to PQ event detection [5]. In the same direction, Gu and Bollen [6] found relevant characteristics associated to PQ events in the time and frequency domains. The work by Ribeiro et al. is also remarkable [7], which extracted new time-domain features based in cumulants. The same authors performed the classification of single and multiple disturbances using HOS in the time domain and Bayes’ theory-based techniques [8]. HOS techniques and estimators have also been implemented to specifically detect sags and swells [9].

The categorization of PQ anomalies had been formerly performed by Nezih and Ece in the work [10], where they proved that HOS and quadratic classifiers improve the second-order-based methods. The same authors previously achieved performance in second-order computing, using 2-D wavelets and compression techniques [11,12]; finding, despite the promising results, the limits of the procedure and quantifying its heavy computational cost. Alienated to this work, the researches of Poisson et al. and Santoso et al. [13,14] also reported a wavelet-based method, finding the potential and the drawbacks of the technique so that to implement it in an intelligent meter.

The direct antecedent of the present research in the work by J.J.G. de la Rosa et al. [15], in which they performed a mixed study involving the time-domain variance, skewness, and kurtosis, and they obtained consequences combined with the spectral kurtosis (SK), and over a set of real-life measurements, some of them with mixed PQ perturbations. The same authors proposed a preliminary criteria for seven types of disturbances based on the former estimators [16]. In a previous work [17], they designed an offline case-based reasoner based on time-domain HOS estimators. Furthermore, the authors used also HOS features in classification techniques for characterization of electrical PQ signals [18].

The present paper is designed as follows. The next section summarizes the main advances in the field, paying special attention at the applications of ANN for PQ analysis, and reasoning the contribution of HOS to the feature extraction stage and the architecture. Then, Section Higher-order statistics for PQ monitoring: an enhancement proposal over the DWT exposes the advantages of HOS for PQ monitoring and the state-of-the art. In Section Proposed methodology: the HOS-based ANN, the procedure is detailed in order to expose the results later in Section Results; finally, conclusions are drawn in Section Conclusions.

2 ANN for PQ analysis: towards the HOS paradigm

The present paper postulates non-stationary signal processing with higher-order statistics in the time and frequency domains in order to extract a battery of features to be processed via ANNs with regression algorithms. ANNs have been used for classification purposes in myriads of works and have proved the utility for a long time [19]; Figure 1 shows a generic architecture, indicating how the input pattern is processed via the neurons to provide with the output vector, in order to introduce the concrete architecture in the present work.

Figure 1

Generic ANN structure which postulates the feature extraction stage previous to the ANN.

The performance of the potential PQ monitoring system and consequently the ANN is directly related to the pre-processing and feature extraction techniques used. The main goal of the feature extraction is to represent the data set in a new feature space in which the probability to distinguish classes is higher than the one in the original space. Therefore, the identification of efficient pre-processing and feature extraction techniques is a key issue [20].

Some works have shown that the best performance for non-stationary signals is achieved combining time and frequency characteristics obtained in the feature extraction stage via wavelet transforms (WT), on the hypotheses that WT performance is suitable for the analysis of transient signals with easy computation. Three works are worth to be cited as they are directly concerned with our research, involving ANN and PQ analysis. In the first one to mention, Angrisani et al. presented a wavelet network for automated transient detection; basically, it was based in an extended perceptron which incorporated wavelet nodes [21]. They proposed a unique structure based on the introduction of wavelet nodes in the traditional ANN, in the neurons of the first layer, and was applied successfully to typical transients in the PQ analysis. The number of nodes was fixed to 12, and this number was related to the time location of the perturbation. The network has to be pre-conditioned in order to detect concrete transient signals.

In the same former line, using the WT as a preprocessing tool for ANN, the second work by Iñigo Monedero et al. [22] developed a neural perceptron-based (three hidden layers) real-time first system prototype, which was trained using a disturbance generator for pseudo-synthetic waveforms. They obtained a performance of 89% success combining simulated and real-time data. The perturbations were classified in frequency, voltage, and harmonics, and the topology of the ANN (number of hidden neurons and number of outputs) was highly moveable depending on the PQ event under study.

The third and most recent paper, by Martin Valtierra-Rodriguez et al. [23] involves a new dual neural-network-based methodology to detect and classify single and combined PQ disturbances, consisting, on the one hand, of an adaptive linear network for harmonic and inter-harmonic estimation that allowed computing the root-mean-square voltage and THD indices, from which it is possible to detect and classify sags, swells, outages, and harmonics. A complementary feed-forward neural network was used for pattern recognition using the horizontal and vertical histograms of a specific voltage waveform to classify spikes, notching, flicker, and oscillatory transients. The complementary action of these neural networks allows the detection and classification even with simultaneous electrical anomalies, both in noisy and noise less scenarios.

Apart from the difficulties in implementing a real-time processor in the smart grid frame, the main disadvantage of the above procedures lies in the fact that accuracy and repeatability are highly compromised by the second-order estimators (WT and V _RMS) used in the feature extraction stage. Drawbacks arise when the data are corrupted by noise; specially when the number of samples of the signal window is reduced, and the resolution and repeatability are degraded. This facts are dramatically increased if the tests are performed over synthetic signals of controlled-lab experiences, where predictability is tacitely supposed. In the real-world experiences, the system should be prepared for unpredictable phenomena, both in the time and in the frequency domains. As explained hereinafter, this goal is accomplished by statistical parameters of an order higher than two.

3 Higher-order statistics for PQ monitoring: an enhancement proposal over the DWT

During the last decade, some researches [7,15,24] have demonstrated the usage of HOS features for PQ monitoring. The motivation of HOS in PQ analysis is twofold. By one side, as HOS measurements are correlations that involve powers higher than two [25], the HOS information for Gaussian signals is null, as Gaussian processes are described up to the second order. As many measurement noises are Gaussian or symmetrically distributed, the HOS may be less degraded by the background noise than the second-order calculations. Secondly, the capability to reveal non-linear characteristics from the data, which is important for pattern recognition scenarios, in concrete to target PQ events, waveforms which exhibit high time variability, like impulsiveness of peakedness.

The usage of HOS as a feature extraction technique for PQ monitoring systems is very promising, and several recent works presented good results with respect to both detection and classification tasks. The most similar work was recently committed by Liu et al. [26]; transients were classified according to the features extracted via the SK and using ANNs. They worked with simulated signals over the five types of synthetic perturbations.

Formerly, it had been shown that combining techniques allows efficient classification of single and simultaneous disturbances, and more, the usage of the second- and fourth-order HOS features, for a specific lag chosen from Fisher’s discriminant ratio (FDR) criterion, has been enough to deal with the majority of the disturbances considered [24].

Regarding specifically signal processing issues, some results show that HOS are capable of detecting disturbances even using short acquisition time windows, which represents an important characteristic for several power system applications such as protection, signal segmentation, and disturbance localization. Being specific, the results shown that the detection of disturbances can be accomplished in less than a quarter of cycle, which is excellent for protection application, where speed and accuracy need to be combined to guarantee selectivity and reliability during the occurrence, for example, of a fault in a system [24,25].

Regarding the HOS that has been used in the present research, it is worth remarking that we have made use of a combination in the time and the frequency domains. The ensemble of Equation 11, whose rigorous mathematical treatment is in the appendix, constitutes indirect measurements of the variance, skewness, and kurtosis that has been used in the present paper, along with the SK, which is used to locate the transients in the frequency domain.

With the aim of motivating the performance of HOS over second-order methods, a comparison has been developed via the discrete wavelet transform (DWT), which is a promising technique for PQ analysis according to the literature (e.g., [21]). Without the further necessity of expanding expressions, in this section, we recall that every finite energy signal s(t) can be decomposed over a wavelet orthogonal basis according to:

$$ s(t) = \sum\limits_{j=-\infty}^{+\infty}\sum\limits_{k=-\infty}^{+\infty}\langle{s,\psi_{j,k}}\rangle\psi_{j,k}. $$

((1))

Each partial sum, indexed by k, in Equation 1 represents the detail variations at the scale a=2^j (at each level j, the scale is increased by a factor of two):

$$ d_{j}(t) = \sum\limits_{k=-\infty}^{+\infty}\langle{s,\psi_{j,k}}\rangle\psi_{j,k} \qquad s(t) = \sum\limits_{j=-\infty}^{+\infty}d_{j}(t). $$

((2))

The approximation of the signal s(t) can be progressively improved by obtaining more levels, with the aim of recovering the signal selectively. For example, if s(t) varies smoothly, we can obtain an approximation by removing fine scale details, which gather information regarding the high frequencies or rapid variations of the signal. This is done by truncating the sum in Equation 1 at the scale a=2^J:

$$ s_{J}(t) = \sum\limits_{j=J}^{+\infty}d_{j}(t). $$

((3))

Details corresponding to indexes j <J are not considered in Equation 3.

Hereinafter, a simple experience is exposed illustrating the limitations of the method based in the DWT, in case DWT coefficients were selected to be used as features. The wavelet decomposition tree has been expanded up to level 7 using the mother wavelet sym8. The idea in the background literature is to study the distribution of the signal energy for each level of the tree, which acts as a filter bank. As stated and proved in myriads of references (e.g., [27]), low index levels point rapid fluctuations of the signal, whereas slow variations are associated to higher (or deeper) decomposition branches. The difference between the energy of the perfect power signal and the energy of the PQ event under test has been studied for each level. The energy expression that we have considered is the traditional 2-norm, according to Equation 4:

$$ E(\textbf{s}) := \sum\limits_{i=1}^{N} \|s_{i}\|^{2}, $$

((4))

where s is the N-point vector signal and s _i , their associated components. Figure 2 shows the processing results for six prototypes. The limitations of the DWT can be concluded via a simple analysis of the graphs. Firstly, as the DWT analysis consists of a selective filtering, sag and swells cannot be distinguished. Secondly, mixed PQ events are also difficult to separate because the constant incremental energy does not indicate whether a transient is coupled or not. Finally, an eventual noisy environment would mask results in the case of the oscillatory transient (this has been extensively proved).

Figure 2

DWT analysis for different PQ events. From top to bottom and left to right: sag, oscillatory transient, sag plus oscillatory transient, swell, harmonics, and pure power signal. The energy difference is depicted vs. the decomposition level. In the case of additive noise, the method does not guarantee success in the feature extraction stage.

At the light of the second-order analysis, a more robust set of features is required in order to process real-life measurements. This is achieved via HOS.

4 Proposed methodology: the HOS-based ANN

The final goal of this work is the implementation in the future SG to guarantee compatibility between all equipment connected. Data used in this work have been simulated using MATLAB software. The data set generated by the simulation consists of 550 samples including the different studied disturbance kinds, which cover the following disturbances: oscillatory transient, impulsive transient, interruption, harmonic permanent distortion, harmonic temporal distortion, sag, sag plus oscillatory transient, and swell. Each signal comprises a 20K-point synthetic time-domain register with a duration of 1 s (20 KHz sampling frequency). An additive normal noise process (1% of the amplitude of the signal) has been added in order to achieve a more realistic behavior. Figure 3 shows an example of these signals as well as healthy signal that are utilized in this paper.

Figure 3

Example of healthy signal and different disturbances studied.

The classification techniques used to classify the PQ disturbances are based on regression algorithms. These techniques are multiple linear regression (MLR), adaptive linear neuron (LIN), multilayer ANNs (BP1 and BP2), radial basis function (RBF) network, exact radial basis (ERB), and generalized regression network (GRN). Every model has been intelligently adapted to meet the objective of PQ classification because each has different characteristics. Table 1 shows the selected parameters corresponding to the architecture and activation function of each model used. The final design of each model has been obtained by the optimization of the parameters shown in this table.

Table 1 Parameters of the ANNs

Data used to realize the classification are based on representative coefficients obtained from the PQ disturbances referred above. These coefficients are acquired by a process of feature extraction which is based on the combination of higher-order statistics in time and frequency domains. The HOS have been computed using a 400-points sliding window (which corresponds to a signal period), with a shift of 10 points over a vector of 20,000 points. After extraction stage, a total of 14 characteristic features are selected, nine of whom correspond to time domain and the remaining five to frequency domain. The coefficients selected in the first one correspond to the maxima and minima and stable in the second-, third-, and fourth-order cumulants at zero lags (directly related to the variance, skewness, and kurtosis). On the other hand, the coefficients selected in the second one correspond to the frequency of extreme value of SK, bandwidth of dome, extreme value of SK, number of peaks in SK, and dome very targeted (between 0 and 1). The abstract graphic of this work with the proposed methodology can be seen in Figure 4 where it specifies the feature extraction module. In order to illustrate the capability of the SK to discriminate PQ disturbances, we have selected a practical example consisting of an oscillatory transient coupled to the power sine wave. The analysis result is depicted in Figure 5. The time-domain variance increases when it bumps into the transient; this behavior is independent of the transient frequency. Similarly, the time-domain skewness and kurtosis detect slight variations. The real detection takes place in the frequency domain; the SK produces a real enhancement in 2,000 Hz, along with the high-resolution bump (narrow peak).

Figure 4

Procedure of feature extraction.

Figure 5

An example of mixed analysis in the time and frequency domains. An oscillatory transient is detected in the time domain by the second-, third-, and fourth-order statistics, and its frequency is clearly targeted by the spectral kurtosis.

The resulting data after carrying out the feature extraction is a matrix of dimension 550×14 (samples × features). Then, the building of the models is performed by following two steps in order to efficiently classify the disturbances. First, data are normalized so that they are in the interval [ −1,1], for a faster computation [28]. And second one, we divide randomly the data set into three subsets: training, evaluation, and test sets. The training and validation sets, with 70% and 15% of data, respectively, were used for ANN model building; and the third set, with the last 15%, was used to test on the out-of-sample set the classification power of a model.

For comparison purposes, the classification criterion is based on the parameter hit rate (HR), that is defined as follows:

$$ \text{{HR}} = \frac{N_{C}}{N_{T}} \times 100\% $$

((5))

where N _T is the number of test samples and N _C is the number of correct disturbance recognition.

5 Results

Once the assessed models were configured and optimized, they are used in the out-of-sample set. Because the database is small, we choose randomly two test sets, and in each of them, 100 experiments are launched by model. This is realized to achieve statistically meaningful results which rule out the random factors influencing the ANNs.

For each model, we have performed three analysis in function of the features used: HOS in the time domain (H O S _t ), the SK, and the mixed (time and frequency domains) analysis (H O S _t +S K). In all of them, the HR was calculated for each model on the 200 offline tests, observing the arithmetic mean. The obtained results showing the percentage of effectiveness for classifying disturbances are presented in Table 2. In most models, results obtained in the mixed analysis are better than in the other analyses except isolated results. These exceptions represent the 27.27% of the cases presented in Table 2, where the 14.28% of them correspond to (H O S _t ) features, and the remaining 12.99% correspond to SK features.

Table 2 Regression algorithms for ANNs

As can be seen, the algorithms based on linear models are worse than those based on non-linear schemes. The best models both individually and collectively are ERB, GRN, and BP1.

6 Conclusions

In this paper, seven regression algorithms have been applied and compared for PQ disturbances classification. The novel aspect is the introduction of new representative coefficients based on HOS in time and frequency domains. These coefficients are the inputs used in the classification algorithms to verify the occurrence or not of single or multiple disturbances in the electric signals.

The data used to test the proposed method were generated by the MATLAB software. The PQ disturbances considered are the most common on the supply. The best models both individually and collectively are obtained employing radial basis networks, generalized regression, and multilayer perceptron. The overall hit rates obtained are 94.70%, 79.59%, and 74.17%, respectively. This is consistent with non-linearity of the used data and emphases of the non-linearity that provide the HOS.

Once the obtained results demonstrate that the proposed method can effectively classify different kinds of PQ disturbances, is necessary to do more training and tests to find an ANN classifier with better characteristics than obtained in this research.

7 Appendix

7.1 Time-domain cumulants and higher-order spectra

7.2 Time-domain HOS

Higher-order cumulants are being used extensively to deduce newly statistical features from the data of non-Gaussian measurement time series [29-33]. To reach a compact expression, nomenclature is first introduced.

Let us consider {x(t)} be an rth-order stationary real-valued random process; the rth-order cumulant is defined as the joint rth-order cumulant of the random variables x(t), x(t+ τ ₁ ),…, x(t+ τ _r−1 ). This compacted notation is expressed via Equation 6:

$$ C_{r,x}(\tau_{1}, \tau_{2},\ldots, \tau_{r-1}) = {\rm{Cum}}[\!x(t),x(t+\tau_{1}),\ldots, x(t+\tau_{r-1})] $$

((6))

where τ ₁,τ ₂,…,τ _r−1 are time shifts, and the nth shifting is a multiple of the data acquisition sampling period, T _s , and is usually expressed as τ _n =n·T _s .

Cumulants, defined in the Equation 6, are estimated by using the well-known Leonov-Shiryaev formula, which expresses the compact relationship among the cumulants of stochastic signals and their moments [9]. In this sense, the expressions for the second-, third-, and fourth-order cumulants for a real, random, and zero-mean (central cumulants) time series x(t) can be estimated via:

$$ C_{2,x}(\tau) = E\{x(t)\cdot x(t + \tau)\} $$

((7a))

$$ C_{3,x}(\tau_{1},\tau_{2}) = E\left\{x(t)\cdot x(t + \tau_{1})\cdot x(t + \tau_{2})\right\} $$

((7b))

$$ \begin{aligned} C_{4,x}(\tau_{1},\tau_{2},\tau_{3})& = E\{x(t)\cdot x(t+ \tau_{1})\cdot x(t+ \tau_{2})\cdot x(t+\tau_{3})\} \\ &\quad- C_{2,x}(\tau_{1})C_{2,x}(\tau_{2}-\tau_{3})\\ &\quad-C_{2,x}(\tau_{2})C_{2,x}(\tau_{3}-\tau_{1})\\ &\quad- C_{2,x}(\tau_{3})C_{2,x}(\tau_{1}-\tau_{2}) \end{aligned} $$

((7c))

where E{∗} is the expected value operator. Then, looking at Equation 7, each cumulant is easily interpreted as a correlation between the original time series and its associated time-shifted versions, being the computational result of an rth-order cumulant is the rth degree of similarity among the aforementioned time series.

Considering a finite N-sample vector, signal vector x(n), n=0,⋯,N−1, the following expressions, Equations 8, 9, and 10 describe the three unbiased estimates for the second-, third-, and fourth-order cumulants, respectively:

$$ \hat{C}_{2,x}(k) = \frac{1}{N} \sum\limits_{n=0}^{N-1} \left[ x(n)\right] \left[x(n+\tau)\right], $$

((8))

$$ \begin{aligned} \hat{C}_{3,x}(k,l) &= \hat{\rm{Cum}}[\!x(n),x(n+k),x(n+l)] \\ &= \frac{1}{N}\sum\limits_{n=0}^{N-1} x(n)x(n+k)x(n+l) \end{aligned} $$

((9))

$$ {\fontsize{8}{6} \begin{aligned} \hat{C}_{4,x}(k,l,m)& = \hat{\rm{Cum}}[\!x(n),x(n+k),x(n+l),x(n+m)]\\ & = \frac{1}{N}\sum\limits_{n=0}^{N-1} x(n)\cdot x(n+k)^{*} \cdot x(n+l)^{*} \cdot x(n+m)^{*}\\ & - \frac{1}{N^{2}}\left[\sum\limits_{n=0}^{N-1} x(n)\cdot x(n+k)^{*} \right] \left[ \sum\limits_{n=0}^{N-1} x(n+l)^{*} \cdot x(n+m)^{*} \right]\\ & - \frac{1}{N^{2}}\left[ \sum\limits_{n=0}^{N-1} x(n)\cdot x(n+l)^{*} \right] \left[ \sum\limits_{n=0}^{N-1} x(n+k)^{*} \cdot x(n+m)^{*} \right]\\ & - \frac{1}{N^{2}}\left[ \sum\limits_{n=0}^{N-1} x(n)\cdot x(n+m)^{*} \right] \left[ \sum\limits_{n=0}^{N-1} x(n+k)^{*} \cdot x(n+l)^{*} \right]\\ \end{aligned}} $$

((10))

where k,l,m∈[ −χ,…,−1,0,1,…,+χ] and n=0,1,…,N−1; χ is the index of the maximum time shift (lag) between samples of a record. The biased expressions are estimates over the real terms in the summations of expressions 8, 9, and 10. These expressions establish the correlation between the original signal and its time-shifted versions for the three orders of comparison. The second-order version is the classical auto-correlation, the third-order one account with the symmetry of the signal, and the fourth-order cumulant quantifies the impulsiveness in the time domain.

Avoided time shifting, τ ₁=τ ₂=τ ₃=0 in Equation 7, leads to the simplest computational expressions for cumulants, in Equation 11:

$$ \gamma_{2,x} = E\{x^{2}(t)\} = C_{2,x}(0) $$

((11a))

$$ \gamma_{3,x} = E\{x^{3}(t)\} = C_{3,x}(0, 0) $$

((11b))

$$ \gamma_{4,x} = E\{x^{4}(t)\}-3(\gamma_{2,x})^{2} = C_{4,x}(0, 0, 0). $$

((11c))

The ensemble of Equation 11 constitutes indirect measurements of the variance, skewness, and kurtosis. If x(t) is symmetrically distributed, its skewness is zero (but not vice versa, improbable situations); if x(t) is Gaussian distributed, its kurtosis is necessarily zero (but not vice versa). Standardization (statistical normalization) makes estimators shift and scale invariant. Standardized quantities are defined as γ _4,x /(γ _2,x )² and γ _3,x /(γ _2,x )^3/2, for kurtosis and skewness, respectively.

7.3 Frequency-domain HOS

Poly-spectra are defined to be the Fourier transforms of the higher-order cumulant sequences. The rth-order spectra are defined as the (r-1)-dimensional Fourier transforms of the rth-order cumulants, according to:

$$ \begin{aligned} S_{r,x}&\left(\,f_{1}, f_{2},\ldots, f_{r-1}\right) =\\ &\sum\limits_{\tau_{1}=-\infty}^{\tau_{1}=+\infty} \dots \sum\limits_{\tau_{r-1}=-\infty}^{\tau_{r-1}=+\infty}C_{r,x}(\tau_{1},\tau_{2},\ldots, \tau_{r-1})\\&\quad\times\exp\left[\!-j2\pi\left(\,f_{1}\tau_{1}+f_{2}\tau_{2}+\dots+f_{r-1}\tau_{r-1}\right)\right]. \end{aligned} $$

((12))

The power spectrum is the decomposition of the signal power in the frequency domain. When this concept is extended to higher orders, as suggested by Equation 12, the result is called a poly-spectrum. Power spectrum, bi-spectrum, and tri-spectrum are specific cases (particular poly-spectra) of Equation 12, with r=2,3, and 4, respectively. Only power spectrum is real, and the others are complex magnitudes.

The more common higher-order spectra are the bi-spectrum and the tri-spectrum. The first one identifies contributions to a signal’s skewness as a function of frequency pairs, meanwhile the tri-spectrum refers to contributions to a signal’s kurtosis as a function of frequency triplets. For this reason, poly-spectra output multidimensional data structures which comprise redundant information, distributed in multi-dimensional geometries, often called tensors. As a consequence, their computation may be impractical in many cases, and to extract the desired information, one-dimensional slices of cumulant sequences and spectra and bi-frequency planes are considered [34,35]. To show this, in the following sections, we present two particular cases of the third- and fourth-order spectra, respectively: a 3-D bi-spectrum application and the performance of an estimator of the fourth-order spectrum for zero time lags, the spectral kurtosis (SK).

Ideally, the spectral kurtosis is a representation of the kurtosis of each frequency component of a process (or data from a measurement instrument x _i ). For estimation issues, we will consider M realizations of the process; each realization containing N points; i.e., we consequently consider M measurement sweeps, each sweep with N points. The time spacing between points is the sampling period, T _s , of the data acquisition unit. The SK unbiased indirect estimator is given by Equation 13:

$$ \widehat{G}_{2,X}^{N,M}=\frac{M}{M-1} \left[\frac{\left(M+1\right)\sum_{i=1}^{M}\left|{X_{N}^{i}}(m)\right|^{4}} {\left(\sum_{i=1}^{M}\left|{X_{N}^{i}}(m)\right|^{2}\right)^{2}}-2\right] $$

((13))

where m indicates the frequency index and $\hat {G}_{2,X}^{N,M}$ indicates the value of the kurtosis for this Fourier frequency. This expression offers an indirect calculation of the SK, as it is obtained directly from the Fourier transforms, and it supposes low computational burden. The graphical representation of the SK allows the identification of non-Gaussian frequency components. The higher the peak the more variable is the amplitude associated to this Fourier component.