- Email: [email protected]
LINEAR OR NONLINEAR? A BICOHERENCE BASED METRIC OF NONLINEARITY MEASURE
LINEAR OR NONLINEAR? A BICOHERENCE BASED METRIC OF NONLINEARITY MEASURE M. A. A. Shoukat Choudhury ∗ David S. Shook ∗ Sirish L. Shah ∗∗,1 ∗
Matrikon Inc., Suite 1800, 10405 Jasper Avenue Edmonton, AB, Canada T5J 3N4 ∗∗ Department of Chemical and Materials Engineering, University of Alberta, Edmonton AB, Canada, T6G 2G6
Abstract: Most process controllers are linear and do perform well when they are optimally tuned. However, the presence of nonlinearities in a process can affect the control performance adversely, e.g., valve stiction may cause oscillations; other nonlinearities e.g. due to process operating conditions change, may render the current tuning parameters as ‘non-optimal’ and so on. For this reason a data-based quantification of process nonlinearity is useful. A new Total Nonlinearity Index (T N LI) is defined to quantify nonlinearities. The method is based on a better estimation of the bicoherence with less number of spurious peaks. The previously published Non-Gaussianity Index (N GI) and Nonlinearity Index (N LI) are modified to overcome their limitations. The efficacy of the proposed indices has been shown using two illustrative examples and one simulation c example. Copyright °2006 IFAC Keywords: Nonlinearity, Bispectrum, higher order statistics, control loops
1. INTRODUCTION
Collis et al. (1998) suggested an addition of a small constant to the denominator to remove these spurious peaks. However, the use of a fixed or stationary constant does not produce reliable results (more discussions at Section 2.2). This work describes a method to choose this constant dynamically depending on the noise level of the time series. Choosing the small constant dynamically helps to obtain a better estimate of the bicoherence with a significant reduction of spurious peaks. The previously published N GI and N LI (Choudhury et al., 2004) have been modified to overcome their limitations. A new total nonlinearity index (T N LI) is introduced as a metric or measure quantify nonlinearities.
Classical signal processing tools utilize only the first and second order moments, i.e., the mean, variance, covariance and correlation. Such tools are mainly useful for analyzing signals from linear processes. The distribution of signals from nonlinear processes is often skewed and non-normal. This necessitates the use of higher order statistical tools. The third and fourth order moments and their frequency domain counterparts (bispectrum and trispectrum) are found to be useful in analyzing nonlinearities in communication signals and mechanical machine condition monitoring (Nikias and Petropulu, 1993; Collis et al., 1998). These higher order statistical techniques have also been used to detect and diagnose nonlinearities in control valves used in process industries (Choudhury et al., 2004; Choudhury, 2004).
2. BISPECTRUM AND BICOHERENCE
In the estimation of the bicoherence, many spurious peaks arise due to the occurrence of small magnitudes of the denominator used to normalize the bispectrum. 1
Corresponding author’s email: [email protected], Tel: +1 (780) 448 1010 ext. 4561, Fax: +1 (780) 448 9191.
617
Industrial processes generally deviate from Gaussianity and linearity and exhibit nonlinear behaviour. These processes can conveniently be studied using Higher Order Statistics (HOS) (Nikias and Petropulu, 1993). Among the various frequency domain HOS measures the bispectrum is the simplest. It is
enough to mask the significant peaks. Therefore, the magnitude of ε needs to be changed dynamically with the calculation of bicoherence for each time series. This point will be illustrated by the following example.
the frequency domain counterpart of the third order moments and defined as B(f1 , f2 ) , E[X(f1 )X(f2 )X ∗ (f1 + f2 )]
(1)
where, X(f ) is the Fourier transformation of the data series x(t). Equation 1 shows that it is a complex quantity having both magnitude and phase. It can be plotted against two independent frequency variables, f1 and f2 in a three dimensional plot. Each point in the plot represents the bispectral energy content of the signal at the bifrequency, (f1 , f2 ). In fact, the bispectrum,B(f1 , f2 ), at point (f1 , f2 ) measures the nonlinear interaction between frequencies f1 and f2 (Nikias and Petropulu, 1993). This interaction between frequencies can be related to the nonlinearities present in the signal generating systems and therein lies the core of its usefulness in the detection and diagnosis of nonlinearities.
2.2 Illustrative Example 1 This example illustrates the effect of ε on the estimation of bicoherence of a Quadratically Phase Coupled (QP C) signal under the influence of various noise levels. Let a QP C signal be constructed as follows: y(t) = sin(2πf1 t + φ1 ) + sin(2πf2 t + φ2 ) + 0.1 sin(2πf3 t + φ3 ) + n(t)
The bispectrum is normalized in the following way to give a measure called bicoherence whose magnitude is bounded between 0 and 1 : bic2 (f1 , f2 ) ,
|E[B(f1 , f2 )]|2 E[|X(f1 )X(f2 )|2 ]E[|X(f1 + f2 )|2 ]
(2)
where ‘bic’ is known as ‘bicoherence’ function. Equation 2 can be rewritten as bic2 (f1 , f2 ) ,
|E[X(f1 )X(f2 )X ∗ (f1 + f2 )]|2 E[|X(f1 )X(f2 )|2 ]E[|X(f1 + f2 )|2 ]
(3)
If Welch’s periodogram method is used to estimate the bicoherence, the expectation operator can be replaced with a summation operator over the number of data segments using the assumption of ergodicity: PM 1 Xi (f1 )Xi (f2 )Xi∗ (f1 + f2 )|2 |M i=1 ˆ 2 (f1 , f2 ) = bic PM P M 1 1 2 M
i=1
|Xi (f1 )Xi (f2 )|
· M
i=1
|Xi (f1 + f2 )|2 (4)
The squared bicoherence is usually estimated using this equation.
2.1 Spurious Peaks in the estimated bicoherence Often times there are small spurious peaks in the bicoherence plot which make the interpretation of the peaks difficult (see the left panel in Figure 2). This Figure is discussed in detail in the next section. Fackrell et al. (1995) reported that these spurious peaks are due to the small values occurring in the denominator of the bicoherence (Equation 4) which cause the estimate to become ill-conditioned. To mitigate such effects, they suggest adding a small constant, ε, to the denominator of Equation 2: |E[B(f1 , f2 )]|2 (5) E[|X(f1 )X(f2 )|2 ]E[|X(f1 + f2 )|2 ] + ε
It is also reported in their paper that the disadvantage of using this correction technique is to introduce a negative bias error in the bicoherence estimator. One needs to select ε such that it is large enough to remove spurious artifacts, but not so large as to cause significant bias effects. However, they did not suggest a suitable numerical value or any method to choose a correct value of ε. In our experience of working with many industrial time series data, a static numerical value of ε does not provide a good estimate of the bicoherence for all cases (see the middle panel in Figure 2). Also, it may introduce bias errors large
where the values of f1 , f2 and f3 are 0.12, 0.18 and 0.30, respectively; the values of φ1 , φ2 , and φ3 are π π 5π 3 , 12 and 12 , respectively; n(t) is a zero-mean white noise signal and t corresponds to discrete sample instants from 1 to 4096 s. The signal y(t) is a quadratic phase coupled signal because its frequencies have the relation f1 + f2 = f3 and its phases have the relation φ1 + φ2 = φ3 . Therefore, the phase coupling at bifrequency (0.12, 0.18) should appear as a single peak in the bicoherence plot. Let the signal to noise ratio (SN R) be defined as the ration of the variance of the noise-free signal to the variance of the noise. All subplots in Figure 2 show that there is one peak at the bifrequency (0.12, 0.18). The top panel of Figure 2 shows the case for SN R = 3 while the bottom panel is for SN R = 10. As evident from this figure, if ε = 0 is used, there are many small spurious peaks in the bicoherence plot. As expected, the scenario is worse for the case of low signal to noise ratio. From the top-left plot of Figure 2, it is hard to find the true peak. For a static value of ε (in this case, chosen as ε = 0.00001), the spurious peaks are absent, but the bias effect on the bicoherence magnitude is very significant. Also, it is difficult to choose a proper static value for ε. In fact, one needs to adopt a ‘trial and error’ procedure. For the known signals whose characteristics are known, a proper value for ε using trial and error technique can be chosen without the fear of loosing the significant peaks. But such is not the case for the practical usage of the bicoherence, where the peaks of bicoherence are unknown and uncertain. Therefore, a better method for choosing ε is sought. The following section suggests a new method of choosing ε dynamically. −5
1
x 10
0.8
Denominator, D
bic2 (f1 , f2 ) ,
(6)
0.6 0.4 0.2 0 0
ε
0.1
0.2
0.3
Frequency
0.4
0.5
Fig. 1. Denominator of Eq. 2 for illustrative example 1 618
ε = 0.00001
εproposed
SNR = 10
SNR = 3
ε=0
Fig. 2. Effect of ε in bicoherence estimation new estimation method with dynamic ε allows us to modify the previously published indices. The modified indices have fewer cases of false positives and better capability of detecting nonlinearities. The following section describes the modification of the indices with a brief overview of the old indices.
2.3 How to choose ε? Let us denote the denominator of Equation 4 as D. This denominator is mainly the product of the power spectrum of the signal. As the power spectrum of most signals consists of few peaks, each column or row of denominator D is also characterized with few peaks. Figure 1 shows the columns of D plotted against frequency f1 for 0 ≤ f1 ≤ 0.5. This figure shows that there are only a few peaks and the rest are small mainly due to noise. The noise level can be denoted as a dotted line in Figure 1. Choosing the magnitude of ε at this level helps to get rid of most of the spurious peaks. In order to obtain the value of ε automatically, it can be chosen as the maximum of the P th percentiles of the columns of D. If it is assumed that there will be a maximum of 25% peaks in each column of D, the value of P can be chosen as 75. The dotted line in Figure 1 is drawn at the maximum of the 75th percentiles of the columns of D. However, depending on the noise level a higher or lower value can be chosen.
A discrete ergodic stationary time series, x(k), is called linear if it can be represented by x(k) =
h(n)e(k − n)
(7)
n=0
where e(k) is a sequence of independent identically distributed random variables with E[e(k)] = 0, σe2 = E[e2 (k)], and µ3 = E[e3 (k)]. For this case, the following frequency domain relationships can be obtained. The power spectrum: P (f ) = σe2 | H(f ) |2 ≡ |X(f )X ∗ (f )|
(8)
and the bispectrum: B(f1 , f2 ) = µ3 H(f1 )H(f2 )H ∗ (f1 + f2 )
The last panel of Figure 2 shows the improvement of the bicoherence estimation using the proposed ε. The use of the proposed dynamic ε could clearly remove most of the spurious peaks without introducing a significant negative bias effect on the bicoherence peak, which can be observed by comparing the maximum peak magnitudes for the cases of ε = 0 (left column plots) and εproposed (right column plots).
where H(f ) =
N −1 X
(9)
h(n)e−if n . Equation 2 can be
n=0
rewritten as bic2 (f1 , f2 ) ,
|E[B(f1 , f2 )]|2 E[|X(f1 )X ∗ (f1 )||X(f2 )X ∗ (f2 )|] · E[|X(f1 + f2 )X ∗ (f1 + f2 )|]
≡
3. TEST OF GAUSSIANITY AND LINEARITY OF A SIGNAL Two indices - Non-Gaussianity Index (N GI) and nonlinearity Index (N LI) - were developed in Choudhury et al. (2004) to test the Gaussianity and linearity of a signal or time series. The indices were developed based on the statistical hypothesis test on the average squared bicoherence in the principal domain (0 < f1 < 0.5, f2 < f1 , and 2f1 + f2 < 1) of the bicoherence. The average bicohernce was used in the formulation of the indices because the estimation of bicoherence was not free from spurious peaks. The
N −1 X
|E[B(f1 , f2 )]|2 E[|P (f1 )||P (f2 )|]E[|P (f1 + f2 )|]
For a linear time series, substituting the expressions from Equations 8 and 9, it can be shown that bic2 (f1 , f2 ) =
619
µ23 σe6
(10)
Equation 10 shows that for any linear signal, x, the squared bicoherence will be independent of the bifrequencies, i.e., a constant in the bifrequency plane. If the squared bicoherence is zero, the signal x is Gaussian because the skewness or µ3 is also zero in such a case. Strictly speaking, such a signal should be called non-skewed with a symmetric distribution
Thus, a signal is non-skewed or Gaussian at a confidence level of α if the N GImodif ied is less than or equal to zero. This index has been defined to facilitate the automation of this decision.
instead of a Gaussian one. However, in this work and also in most of the HOS literature (Nikias and Petropulu, 1993; Hinich, 1982; Rao and Gabr, 1980; Collis et al., 1998) the two terms – nonskewed and Gaussian – have been used interchangeably. To check whether the squared bicoherence is constant, two tests are required. One is to determine whether the squared bicoherence is zero, which would show that the signal is Gaussian and thereby the signal generating process is linear. The other is to test for a non-zero constant value of the squared bicoherence which would show that the signal is non-Gaussian but the signal generating process is linear.
3.2 Modified Nonlinearity Index (N LImodif ied ) If a signal is found to be Gaussian, the signal generating process is assumed to be linear (Rao and Gabr, 1980). In the case of a non-Gaussian signal, the signal generating process should be tested for its linearity. As shown in Equation 10, if the signal is non-Gaussian and linear, the magnitude of the squared bicoherence should be a non-zero constant at all bifrequencies in the principal domain because µ is a non-zero constant.
3.1 Modified Non-Gaussianity Index (N GImodif ied ) The bicoherence is a complex quantity with real and imaginary parts. The magnitude of the squared bicoherence can be obtained as bic2 = <(bic)2 + =(bic)2
A simple way to confirm the constancy of squared bicoherence is to have a look at the 3-D squared bicoherence plot and observe the flatness of the plot. If the squared bicoherence is of a constant magnitude at all bifrequencies in the principal domain, the variance of the estimated bicoherence should be zero. To check the flatness of the plot or the constancy of the squared bicoherence, Choudhury et al. (2004) developed a nonlinearity index by comparing the maximum squared bicoherence with the average squared bicoherence plus two times the standard deviation of the estimated squared bicoherence. The disadvantage of using this index is that the presence of a few large peaks significantly bias the standard deviation and mean of the estimated bicoherence, which leads to some false negatives. The other problem with this index is that it cannot give an indication of the extent of nonlinearity because with the increase of peak sizes, the mean and the standard deviation also increase. In order to avoid these limitations, the index can be modified as:
(11)
where < and = are real and imaginary parts, respectively. It is well established in the HOS literature that bicoherence is a complex normal variable, i.e., both estimates of real and imaginary parts of the bicoherence are normally distributed (Hinich, 1982) and asymptotically independent, i.e., the estimate at a particular bifrequency is independent of the estimates of its neighboring bifrequencies. Therefore, the squared bicoherence at each bifrequency is a non-central chisquared (χ2 (m)) distributed variable with 2 degrees of freedom and mean, m. The following statistical test can be applied to check for the significance of bicoherence magnitude at each individual bifrequency: 2
P {2K bic2 (f1 , f2 ) > cχα } = α or,
ˆ 2 max − (bic ˆ2 N LImodif ied , bic robust + 2σbic ˆ 2 , robust ) (15) 2
ˆ where bic ˆ 2 , robust are the robust robust and σbic mean and the robust standard deviation of the estimated squared bicoherence. They are calculated by excluding the largest and smallest Q% of the bicoherence. A good value of Q may be chosen as 10. There is a similar concept used in statistics to calculate robust mean called ‘trimmed mean’. The users may choose a different value of Q.
(12)
2
P {bic2 (f1 , f2 ) >
cχα }=α 2K
(13)
where K is the number of data segments used 2 in bicoherence estimation, cχα is the critical value 2 calculated from the central χ distribution table for a significance level of α at degrees of freedom 2. For 2 example, at α = 0.05, the value of cχ0.05 is 5.99.
Therefore, it can be concluded that: • if N LImodif ied ≤ 0, the signal generating process is LINEAR • if N LImodif ied > 0, the signal generating process is NONLINEAR
Often the principal domain of the bicoherence plot contains more than a hundred bifrequencies. The hypothesis test results for this large number of bifrequencies can be conveniently summarized into the following modified Non-Gaussianity Index (N GImodif ied ): P 2 2 N GImodif ied ,
bicsignif icant L
−
cχ α 2K
Since the squared bicoherence is bounded between 0 and 1, the Nonlinearity Index (N LI) is also bounded between -1 and 1.
(14)
where bic2signif icant are those bicoherence which fail the hypothesis test in Equation (13), i.e., 2
3.3 Total Nonlinearity Index (T N LI)
2
cχ α 2K
, and L is the number of bic (f1 , f2 ) > bic2signif icant . Therefore, the following rule-based decision is suggested: • if N GImodif ied ≤ α, the signal is GAUSSIAN • if N GImodif ied > α, the signal is NONGAUSSIAN
620
It is important to measure nonlinearity in terms of a metric especially when there is a need to compare the extent of nonlinearities in various time series. There are only a few time domain methods available in the literature to measure the nonlinearity. Thornhill et al. (2003) and Helbig et al. (2000) measure the
nonlinearity in terms of a metric in time domain. However, to the best knowledge of the authors, there is no such metric in the frequency domain. This work proposes a new metric to measure nonlinearity using frequency domain bicoherence. If a time series is detected as nonlinear by the above-mentioned N GImodif ied and N LImodif ied indices, then the total nonlinearity present in the time series can be quantified using the following new index: X T N LI , bic2signif icant (16)
Total Nonlinearity Index, TNLI
6
is Total Nonlinearity Index, bic2signif icant
where, T N LI are those bicoherence which fail the hypothesis test 2 in Equation (12), i.e., 2K bic2 (f1 , f2 ) > cχα . The T N LI is bounded between 0 and L, where L is the number of bic2signif icant .
x (k) = sin(2πf1 k + φ1 ) + sin(2πf2 k + φ2 ) 0
SNR=2 SNR=5 SNR=10
1 0 0
0.2
(17)
where, f1 = 0.12, f2 = 0.30, φ1 = π/3, φ2 = π/8, nl is a multiplication factor employed to represent the contribution of the nonlinear component of the signal, and d(k) is a white noise sequence. Again, frequencies are normalized such that the sampling frequency is 1. The quadratic term in Equation 17 will introduce phase coupled nonlinearity in the output signal, y. The quadratically phase coupled nonlinearities arise due to the interactions between any two of the signal components at frequencies f1 , f2 , 2f1 , 2f2 , f2 − f1 , and f1 + f2 . For details refer to (Choudhury, 2004).
nl
0.6
0.8
1
Fi
R
1 Modified NGI and NLI
0.4
In linear system theory, it is often assumed that a nonlinear process operates in a locally linear fashion. This section investigates the validity of this assumption using the newly defined ‘total nonlinearity index’ through an amplitude and frequency dependent nonlinearity analysis. This analysis is performed using a simulation of the level control of a spherical tank (Agrawal and Lakshminarayanan, 2003). Consider the control of water level in a spherical tank by manipulating the inlet volumetric flow rate, Fi (see Figure 5).
x(k) = x (k) + d(k) y(k) = x (k) + nl x (k) + d(k)
2
5. QUANTIFYING NONLINEARITY OF A SPHERICAL TANK
0
2
3
in Figure 4 represents the extent of nonlinearity in the signal. For nl = 0, there is no nonlinearity. Therefore, T N LI = 0. With increase of nl , the T N LI also increases. The amount of noise influences the calculation of all indices because the denominator of the bicoherence is affected by the noise.
The purpose of this example is to demonstrate the efficacy of the proposed indices in the presence of varying noise and extent of nonlinearity of a signal. Let a signal be generated as follows:
0
4
Fig. 4. T N LI increases with the increase of the nonlinearity of the signal.
4. ILLUSTRATIVE EXAMPLE 2
0
5
0.8
h ho
0.6 0.4
Fig. 5. A spherical tank system. SNR=2 SNR=5 SNR=10
0.2 0 0
Fo
0.2
0.4
nl
0.6
0.8
The dynamics of the system can be modelled as: 1
Fig. 3. N GImodif ied (dotted lines) and N LImodif ied (solid lines). Figure 3 shows the modified Non-Gaussianity Index (N GImodif ied ) and the modified Nonlinearity Index (N LImodif ied ) plotted against nl for varying cases of signal to noise (SN R) ratio. The dotted lines represent the modified N GI and the solid lines show the modified N LI. It is clear from the figure that both indices increase with the increase of nonlinearity in the signal. The total nonlinearity index shown
(R − h)2 dh = Fi (t − d) − Fo (t) (18) ] dt R2 where Fo (t) is the outlet flow rate at time t, h is the height of the water level from the bottom of the tank, R is the radius of the spherical tank, and d is the time delay between Fi and Fo . The poutlet flow rate, Fo , can be expressed as Fo (t) = 2g(h − ho ) where g is the gravitational constant and ho is the height of the outlet pipe from the base of the vessel. For both open loop and closed loop simulations of the system, we use R = 0.5 m, ho = 0.01 m, and the nominal operating point is selected to be, hss = 0.40 m. The steadystate operating point for the water height of the tank πR2 [1 −
621
to investigate this, this process was excited with a sinusoid of amplitude 0.1 m and 0.05 cycle/sample or 20 samples/cycle signal for various operating point. The results are shown in Figure 7. It is clear from the figure that nonlinearity increases as the operating point moves from the center of the tank, i.e., 0.5 m. This is also expected because the curvature of the inner surface of the tank increases as it moves away from the center.
is chosen at 0.4m instead of 0.5m because the process is less nonlinear in the vicinity of the latter operating point. 5.1 Dependency of nonlinearity on the amplitude and frequency of input signal All linear system exhibits ‘sinusoidal fidelity’. Sinusoidal fidelity states that a sinusoidal input to a linear system will produce a sinusoidal output of identical frequency. The amplitude and phase of the output sinusoid may differ from the input sinusoid. But this is not the case for a nonlinear system. Generally, in response to a sinusoidal excitation a nonlinear system will produce new frequencies in addition to the original frequency of the input sinusoid. This property requires the use of a sinusoid as an excitation signal for testing nonlinearity of a system. In this example, a sinusoidal excitation (added as a deviation to the steadystate Fi ) with varying frequencies and magnitudes are used. The amplitudes vary from 0.01 to 0.2 and the frequencies vary from 0.015 to 0.1 cycles/sample or 54 smaples/cycle to 10 samples/cycle. 4096 data points were used for each nonlinear analysis. The results are shown in Figure 6. This figure shows that for small magnitude input signals, the process can be assumed locally linear. As the amplitude of excitation goes up, so does the nonlinearity. Change of frequency of the input signal is not very sensitive to the nonlinearity analysis. Therefore, it can be concluded that nonlinearity of the process strongly depends on the size of the excitation or input signal.
6. CONCLUSIONS A method for better estimation of bicoherence has been presented. The modified N GI and N LI perform better in detecting signal non-Gaussianity and nonlinearity. A new index, Total Nonlinearity Index (T N LI), has been developed and its performance in quantifying signal nonlinearity is demonstrated using illustrative simulation examples. These indices are applicable to any time series signal and does not require a sinusoidal excitation of a system. Further research is in progress towards the use of these bicoherence based indices for root cause diagnosis of plant-wide oscillations. REFERENCES
TNLI
4 3 2 1 0 0.1
0.2 0.08
0.06 0.04 w, frequency 0.02
0.15 0.05
0.1 A, amplitude
0 0
Fig. 6. Dependency of nonlinearity on the input signal’s amplitude and frequency 3.5 3
TNLI
2.5 2 1.5 1 0.5 0
0.2
0.4 0.6 Nominal Level
0.8
1
Fig. 7. Dependency of nonlinearity on the operating region of the system 5.2 Operating point Dependency of nonlinearity It is well-known that the nonlinearity of a system depends on the operating point of the system. In order
622
Agrawal, P. and S. Lakshminarayanan (2003). Tuning pid controllers using achievable performance indices. Ind. Eng. Chem. Res. 42, 5576–5582. Choudhury, M. A. A. S. (2004). Detection and diagnosis of control loop nonlinearities, valve stiction and data compression. PhD thesis. University of Alberta, Canada. Choudhury, M. A. A. S., S. L. Shah and N. F. Thornhill (2004). Diagnosis of poor control loop performance using higher order statistics. Automatica 40(10), 1719–1728. Collis, W.B., P.R. White and J.K. Hammond (1998). Higher-order spectra: The bispectrum and trispectrum. Mechanical Systems and Signal Processing 12, 375–394. Fackrell, J.W.A., P.R. White, J.K. Hammond, R.J. Pinnigton and A. T. Parsons (1995). The interpretation of the bispectra of vibration signals - i. theory. Mechanical Systems and Signal Processing 9, 257–266. Helbig, A., W. Marquardt and F. Allg¨ower (2000). Nonlinearity measure: definition, computation and applications. Journal of Process Control 10, 113123. Hinich, M. J. (1982). Testing for gaussianity and linearity of a stationary time series. Journal of time series analysis 3, 169–176. Nikias, C. L. and A. P. Petropulu (1993). HigherOrder Spectra: A nonlinear signal processing frame work. Prentice Hall. New Jeresy. Rao, T. S. and M. M. Gabr (1980). A test for linearity and stationarity of time series. Journal of time series analysis 1(1), 145–158. Thornhill, N. F., F. Cox and M. Paulonis (2003). Diagnosis of plant-wide oscillation through data-driven analysis and process understanding. Control Engng. Prac. 11(12), 1481–1490.
Matrikon Inc., Suite 1800, 10405 Jasper Avenue Edmonton, AB, Canada T5J 3N4 ∗∗ Department of Chemical and Materials Engineering, University of Alberta, Edmonton AB, Canada, T6G 2G6
Abstract: Most process controllers are linear and do perform well when they are optimally tuned. However, the presence of nonlinearities in a process can affect the control performance adversely, e.g., valve stiction may cause oscillations; other nonlinearities e.g. due to process operating conditions change, may render the current tuning parameters as ‘non-optimal’ and so on. For this reason a data-based quantification of process nonlinearity is useful. A new Total Nonlinearity Index (T N LI) is defined to quantify nonlinearities. The method is based on a better estimation of the bicoherence with less number of spurious peaks. The previously published Non-Gaussianity Index (N GI) and Nonlinearity Index (N LI) are modified to overcome their limitations. The efficacy of the proposed indices has been shown using two illustrative examples and one simulation c example. Copyright °2006 IFAC Keywords: Nonlinearity, Bispectrum, higher order statistics, control loops
1. INTRODUCTION
Collis et al. (1998) suggested an addition of a small constant to the denominator to remove these spurious peaks. However, the use of a fixed or stationary constant does not produce reliable results (more discussions at Section 2.2). This work describes a method to choose this constant dynamically depending on the noise level of the time series. Choosing the small constant dynamically helps to obtain a better estimate of the bicoherence with a significant reduction of spurious peaks. The previously published N GI and N LI (Choudhury et al., 2004) have been modified to overcome their limitations. A new total nonlinearity index (T N LI) is introduced as a metric or measure quantify nonlinearities.
Classical signal processing tools utilize only the first and second order moments, i.e., the mean, variance, covariance and correlation. Such tools are mainly useful for analyzing signals from linear processes. The distribution of signals from nonlinear processes is often skewed and non-normal. This necessitates the use of higher order statistical tools. The third and fourth order moments and their frequency domain counterparts (bispectrum and trispectrum) are found to be useful in analyzing nonlinearities in communication signals and mechanical machine condition monitoring (Nikias and Petropulu, 1993; Collis et al., 1998). These higher order statistical techniques have also been used to detect and diagnose nonlinearities in control valves used in process industries (Choudhury et al., 2004; Choudhury, 2004).
2. BISPECTRUM AND BICOHERENCE
In the estimation of the bicoherence, many spurious peaks arise due to the occurrence of small magnitudes of the denominator used to normalize the bispectrum. 1
Corresponding author’s email: [email protected], Tel: +1 (780) 448 1010 ext. 4561, Fax: +1 (780) 448 9191.
617
Industrial processes generally deviate from Gaussianity and linearity and exhibit nonlinear behaviour. These processes can conveniently be studied using Higher Order Statistics (HOS) (Nikias and Petropulu, 1993). Among the various frequency domain HOS measures the bispectrum is the simplest. It is
enough to mask the significant peaks. Therefore, the magnitude of ε needs to be changed dynamically with the calculation of bicoherence for each time series. This point will be illustrated by the following example.
the frequency domain counterpart of the third order moments and defined as B(f1 , f2 ) , E[X(f1 )X(f2 )X ∗ (f1 + f2 )]
(1)
where, X(f ) is the Fourier transformation of the data series x(t). Equation 1 shows that it is a complex quantity having both magnitude and phase. It can be plotted against two independent frequency variables, f1 and f2 in a three dimensional plot. Each point in the plot represents the bispectral energy content of the signal at the bifrequency, (f1 , f2 ). In fact, the bispectrum,B(f1 , f2 ), at point (f1 , f2 ) measures the nonlinear interaction between frequencies f1 and f2 (Nikias and Petropulu, 1993). This interaction between frequencies can be related to the nonlinearities present in the signal generating systems and therein lies the core of its usefulness in the detection and diagnosis of nonlinearities.
2.2 Illustrative Example 1 This example illustrates the effect of ε on the estimation of bicoherence of a Quadratically Phase Coupled (QP C) signal under the influence of various noise levels. Let a QP C signal be constructed as follows: y(t) = sin(2πf1 t + φ1 ) + sin(2πf2 t + φ2 ) + 0.1 sin(2πf3 t + φ3 ) + n(t)
The bispectrum is normalized in the following way to give a measure called bicoherence whose magnitude is bounded between 0 and 1 : bic2 (f1 , f2 ) ,
|E[B(f1 , f2 )]|2 E[|X(f1 )X(f2 )|2 ]E[|X(f1 + f2 )|2 ]
(2)
where ‘bic’ is known as ‘bicoherence’ function. Equation 2 can be rewritten as bic2 (f1 , f2 ) ,
|E[X(f1 )X(f2 )X ∗ (f1 + f2 )]|2 E[|X(f1 )X(f2 )|2 ]E[|X(f1 + f2 )|2 ]
(3)
If Welch’s periodogram method is used to estimate the bicoherence, the expectation operator can be replaced with a summation operator over the number of data segments using the assumption of ergodicity: PM 1 Xi (f1 )Xi (f2 )Xi∗ (f1 + f2 )|2 |M i=1 ˆ 2 (f1 , f2 ) = bic PM P M 1 1 2 M
i=1
|Xi (f1 )Xi (f2 )|
· M
i=1
|Xi (f1 + f2 )|2 (4)
The squared bicoherence is usually estimated using this equation.
2.1 Spurious Peaks in the estimated bicoherence Often times there are small spurious peaks in the bicoherence plot which make the interpretation of the peaks difficult (see the left panel in Figure 2). This Figure is discussed in detail in the next section. Fackrell et al. (1995) reported that these spurious peaks are due to the small values occurring in the denominator of the bicoherence (Equation 4) which cause the estimate to become ill-conditioned. To mitigate such effects, they suggest adding a small constant, ε, to the denominator of Equation 2: |E[B(f1 , f2 )]|2 (5) E[|X(f1 )X(f2 )|2 ]E[|X(f1 + f2 )|2 ] + ε
It is also reported in their paper that the disadvantage of using this correction technique is to introduce a negative bias error in the bicoherence estimator. One needs to select ε such that it is large enough to remove spurious artifacts, but not so large as to cause significant bias effects. However, they did not suggest a suitable numerical value or any method to choose a correct value of ε. In our experience of working with many industrial time series data, a static numerical value of ε does not provide a good estimate of the bicoherence for all cases (see the middle panel in Figure 2). Also, it may introduce bias errors large
where the values of f1 , f2 and f3 are 0.12, 0.18 and 0.30, respectively; the values of φ1 , φ2 , and φ3 are π π 5π 3 , 12 and 12 , respectively; n(t) is a zero-mean white noise signal and t corresponds to discrete sample instants from 1 to 4096 s. The signal y(t) is a quadratic phase coupled signal because its frequencies have the relation f1 + f2 = f3 and its phases have the relation φ1 + φ2 = φ3 . Therefore, the phase coupling at bifrequency (0.12, 0.18) should appear as a single peak in the bicoherence plot. Let the signal to noise ratio (SN R) be defined as the ration of the variance of the noise-free signal to the variance of the noise. All subplots in Figure 2 show that there is one peak at the bifrequency (0.12, 0.18). The top panel of Figure 2 shows the case for SN R = 3 while the bottom panel is for SN R = 10. As evident from this figure, if ε = 0 is used, there are many small spurious peaks in the bicoherence plot. As expected, the scenario is worse for the case of low signal to noise ratio. From the top-left plot of Figure 2, it is hard to find the true peak. For a static value of ε (in this case, chosen as ε = 0.00001), the spurious peaks are absent, but the bias effect on the bicoherence magnitude is very significant. Also, it is difficult to choose a proper static value for ε. In fact, one needs to adopt a ‘trial and error’ procedure. For the known signals whose characteristics are known, a proper value for ε using trial and error technique can be chosen without the fear of loosing the significant peaks. But such is not the case for the practical usage of the bicoherence, where the peaks of bicoherence are unknown and uncertain. Therefore, a better method for choosing ε is sought. The following section suggests a new method of choosing ε dynamically. −5
1
x 10
0.8
Denominator, D
bic2 (f1 , f2 ) ,
(6)
0.6 0.4 0.2 0 0
ε
0.1
0.2
0.3
Frequency
0.4
0.5
Fig. 1. Denominator of Eq. 2 for illustrative example 1 618
ε = 0.00001
εproposed
SNR = 10
SNR = 3
ε=0
Fig. 2. Effect of ε in bicoherence estimation new estimation method with dynamic ε allows us to modify the previously published indices. The modified indices have fewer cases of false positives and better capability of detecting nonlinearities. The following section describes the modification of the indices with a brief overview of the old indices.
2.3 How to choose ε? Let us denote the denominator of Equation 4 as D. This denominator is mainly the product of the power spectrum of the signal. As the power spectrum of most signals consists of few peaks, each column or row of denominator D is also characterized with few peaks. Figure 1 shows the columns of D plotted against frequency f1 for 0 ≤ f1 ≤ 0.5. This figure shows that there are only a few peaks and the rest are small mainly due to noise. The noise level can be denoted as a dotted line in Figure 1. Choosing the magnitude of ε at this level helps to get rid of most of the spurious peaks. In order to obtain the value of ε automatically, it can be chosen as the maximum of the P th percentiles of the columns of D. If it is assumed that there will be a maximum of 25% peaks in each column of D, the value of P can be chosen as 75. The dotted line in Figure 1 is drawn at the maximum of the 75th percentiles of the columns of D. However, depending on the noise level a higher or lower value can be chosen.
A discrete ergodic stationary time series, x(k), is called linear if it can be represented by x(k) =
h(n)e(k − n)
(7)
n=0
where e(k) is a sequence of independent identically distributed random variables with E[e(k)] = 0, σe2 = E[e2 (k)], and µ3 = E[e3 (k)]. For this case, the following frequency domain relationships can be obtained. The power spectrum: P (f ) = σe2 | H(f ) |2 ≡ |X(f )X ∗ (f )|
(8)
and the bispectrum: B(f1 , f2 ) = µ3 H(f1 )H(f2 )H ∗ (f1 + f2 )
The last panel of Figure 2 shows the improvement of the bicoherence estimation using the proposed ε. The use of the proposed dynamic ε could clearly remove most of the spurious peaks without introducing a significant negative bias effect on the bicoherence peak, which can be observed by comparing the maximum peak magnitudes for the cases of ε = 0 (left column plots) and εproposed (right column plots).
where H(f ) =
N −1 X
(9)
h(n)e−if n . Equation 2 can be
n=0
rewritten as bic2 (f1 , f2 ) ,
|E[B(f1 , f2 )]|2 E[|X(f1 )X ∗ (f1 )||X(f2 )X ∗ (f2 )|] · E[|X(f1 + f2 )X ∗ (f1 + f2 )|]
≡
3. TEST OF GAUSSIANITY AND LINEARITY OF A SIGNAL Two indices - Non-Gaussianity Index (N GI) and nonlinearity Index (N LI) - were developed in Choudhury et al. (2004) to test the Gaussianity and linearity of a signal or time series. The indices were developed based on the statistical hypothesis test on the average squared bicoherence in the principal domain (0 < f1 < 0.5, f2 < f1 , and 2f1 + f2 < 1) of the bicoherence. The average bicohernce was used in the formulation of the indices because the estimation of bicoherence was not free from spurious peaks. The
N −1 X
|E[B(f1 , f2 )]|2 E[|P (f1 )||P (f2 )|]E[|P (f1 + f2 )|]
For a linear time series, substituting the expressions from Equations 8 and 9, it can be shown that bic2 (f1 , f2 ) =
619
µ23 σe6
(10)
Equation 10 shows that for any linear signal, x, the squared bicoherence will be independent of the bifrequencies, i.e., a constant in the bifrequency plane. If the squared bicoherence is zero, the signal x is Gaussian because the skewness or µ3 is also zero in such a case. Strictly speaking, such a signal should be called non-skewed with a symmetric distribution
Thus, a signal is non-skewed or Gaussian at a confidence level of α if the N GImodif ied is less than or equal to zero. This index has been defined to facilitate the automation of this decision.
instead of a Gaussian one. However, in this work and also in most of the HOS literature (Nikias and Petropulu, 1993; Hinich, 1982; Rao and Gabr, 1980; Collis et al., 1998) the two terms – nonskewed and Gaussian – have been used interchangeably. To check whether the squared bicoherence is constant, two tests are required. One is to determine whether the squared bicoherence is zero, which would show that the signal is Gaussian and thereby the signal generating process is linear. The other is to test for a non-zero constant value of the squared bicoherence which would show that the signal is non-Gaussian but the signal generating process is linear.
3.2 Modified Nonlinearity Index (N LImodif ied ) If a signal is found to be Gaussian, the signal generating process is assumed to be linear (Rao and Gabr, 1980). In the case of a non-Gaussian signal, the signal generating process should be tested for its linearity. As shown in Equation 10, if the signal is non-Gaussian and linear, the magnitude of the squared bicoherence should be a non-zero constant at all bifrequencies in the principal domain because µ is a non-zero constant.
3.1 Modified Non-Gaussianity Index (N GImodif ied ) The bicoherence is a complex quantity with real and imaginary parts. The magnitude of the squared bicoherence can be obtained as bic2 = <(bic)2 + =(bic)2
A simple way to confirm the constancy of squared bicoherence is to have a look at the 3-D squared bicoherence plot and observe the flatness of the plot. If the squared bicoherence is of a constant magnitude at all bifrequencies in the principal domain, the variance of the estimated bicoherence should be zero. To check the flatness of the plot or the constancy of the squared bicoherence, Choudhury et al. (2004) developed a nonlinearity index by comparing the maximum squared bicoherence with the average squared bicoherence plus two times the standard deviation of the estimated squared bicoherence. The disadvantage of using this index is that the presence of a few large peaks significantly bias the standard deviation and mean of the estimated bicoherence, which leads to some false negatives. The other problem with this index is that it cannot give an indication of the extent of nonlinearity because with the increase of peak sizes, the mean and the standard deviation also increase. In order to avoid these limitations, the index can be modified as:
(11)
where < and = are real and imaginary parts, respectively. It is well established in the HOS literature that bicoherence is a complex normal variable, i.e., both estimates of real and imaginary parts of the bicoherence are normally distributed (Hinich, 1982) and asymptotically independent, i.e., the estimate at a particular bifrequency is independent of the estimates of its neighboring bifrequencies. Therefore, the squared bicoherence at each bifrequency is a non-central chisquared (χ2 (m)) distributed variable with 2 degrees of freedom and mean, m. The following statistical test can be applied to check for the significance of bicoherence magnitude at each individual bifrequency: 2
P {2K bic2 (f1 , f2 ) > cχα } = α or,
ˆ 2 max − (bic ˆ2 N LImodif ied , bic robust + 2σbic ˆ 2 , robust ) (15) 2
ˆ where bic ˆ 2 , robust are the robust robust and σbic mean and the robust standard deviation of the estimated squared bicoherence. They are calculated by excluding the largest and smallest Q% of the bicoherence. A good value of Q may be chosen as 10. There is a similar concept used in statistics to calculate robust mean called ‘trimmed mean’. The users may choose a different value of Q.
(12)
2
P {bic2 (f1 , f2 ) >
cχα }=α 2K
(13)
where K is the number of data segments used 2 in bicoherence estimation, cχα is the critical value 2 calculated from the central χ distribution table for a significance level of α at degrees of freedom 2. For 2 example, at α = 0.05, the value of cχ0.05 is 5.99.
Therefore, it can be concluded that: • if N LImodif ied ≤ 0, the signal generating process is LINEAR • if N LImodif ied > 0, the signal generating process is NONLINEAR
Often the principal domain of the bicoherence plot contains more than a hundred bifrequencies. The hypothesis test results for this large number of bifrequencies can be conveniently summarized into the following modified Non-Gaussianity Index (N GImodif ied ): P 2 2 N GImodif ied ,
bicsignif icant L
−
cχ α 2K
Since the squared bicoherence is bounded between 0 and 1, the Nonlinearity Index (N LI) is also bounded between -1 and 1.
(14)
where bic2signif icant are those bicoherence which fail the hypothesis test in Equation (13), i.e., 2
3.3 Total Nonlinearity Index (T N LI)
2
cχ α 2K
, and L is the number of bic (f1 , f2 ) > bic2signif icant . Therefore, the following rule-based decision is suggested: • if N GImodif ied ≤ α, the signal is GAUSSIAN • if N GImodif ied > α, the signal is NONGAUSSIAN
620
It is important to measure nonlinearity in terms of a metric especially when there is a need to compare the extent of nonlinearities in various time series. There are only a few time domain methods available in the literature to measure the nonlinearity. Thornhill et al. (2003) and Helbig et al. (2000) measure the
nonlinearity in terms of a metric in time domain. However, to the best knowledge of the authors, there is no such metric in the frequency domain. This work proposes a new metric to measure nonlinearity using frequency domain bicoherence. If a time series is detected as nonlinear by the above-mentioned N GImodif ied and N LImodif ied indices, then the total nonlinearity present in the time series can be quantified using the following new index: X T N LI , bic2signif icant (16)
Total Nonlinearity Index, TNLI
6
is Total Nonlinearity Index, bic2signif icant
where, T N LI are those bicoherence which fail the hypothesis test 2 in Equation (12), i.e., 2K bic2 (f1 , f2 ) > cχα . The T N LI is bounded between 0 and L, where L is the number of bic2signif icant .
x (k) = sin(2πf1 k + φ1 ) + sin(2πf2 k + φ2 ) 0
SNR=2 SNR=5 SNR=10
1 0 0
0.2
(17)
where, f1 = 0.12, f2 = 0.30, φ1 = π/3, φ2 = π/8, nl is a multiplication factor employed to represent the contribution of the nonlinear component of the signal, and d(k) is a white noise sequence. Again, frequencies are normalized such that the sampling frequency is 1. The quadratic term in Equation 17 will introduce phase coupled nonlinearity in the output signal, y. The quadratically phase coupled nonlinearities arise due to the interactions between any two of the signal components at frequencies f1 , f2 , 2f1 , 2f2 , f2 − f1 , and f1 + f2 . For details refer to (Choudhury, 2004).
nl
0.6
0.8
1
Fi
R
1 Modified NGI and NLI
0.4
In linear system theory, it is often assumed that a nonlinear process operates in a locally linear fashion. This section investigates the validity of this assumption using the newly defined ‘total nonlinearity index’ through an amplitude and frequency dependent nonlinearity analysis. This analysis is performed using a simulation of the level control of a spherical tank (Agrawal and Lakshminarayanan, 2003). Consider the control of water level in a spherical tank by manipulating the inlet volumetric flow rate, Fi (see Figure 5).
x(k) = x (k) + d(k) y(k) = x (k) + nl x (k) + d(k)
2
5. QUANTIFYING NONLINEARITY OF A SPHERICAL TANK
0
2
3
in Figure 4 represents the extent of nonlinearity in the signal. For nl = 0, there is no nonlinearity. Therefore, T N LI = 0. With increase of nl , the T N LI also increases. The amount of noise influences the calculation of all indices because the denominator of the bicoherence is affected by the noise.
The purpose of this example is to demonstrate the efficacy of the proposed indices in the presence of varying noise and extent of nonlinearity of a signal. Let a signal be generated as follows:
0
4
Fig. 4. T N LI increases with the increase of the nonlinearity of the signal.
4. ILLUSTRATIVE EXAMPLE 2
0
5
0.8
h ho
0.6 0.4
Fig. 5. A spherical tank system. SNR=2 SNR=5 SNR=10
0.2 0 0
Fo
0.2
0.4
nl
0.6
0.8
The dynamics of the system can be modelled as: 1
Fig. 3. N GImodif ied (dotted lines) and N LImodif ied (solid lines). Figure 3 shows the modified Non-Gaussianity Index (N GImodif ied ) and the modified Nonlinearity Index (N LImodif ied ) plotted against nl for varying cases of signal to noise (SN R) ratio. The dotted lines represent the modified N GI and the solid lines show the modified N LI. It is clear from the figure that both indices increase with the increase of nonlinearity in the signal. The total nonlinearity index shown
(R − h)2 dh = Fi (t − d) − Fo (t) (18) ] dt R2 where Fo (t) is the outlet flow rate at time t, h is the height of the water level from the bottom of the tank, R is the radius of the spherical tank, and d is the time delay between Fi and Fo . The poutlet flow rate, Fo , can be expressed as Fo (t) = 2g(h − ho ) where g is the gravitational constant and ho is the height of the outlet pipe from the base of the vessel. For both open loop and closed loop simulations of the system, we use R = 0.5 m, ho = 0.01 m, and the nominal operating point is selected to be, hss = 0.40 m. The steadystate operating point for the water height of the tank πR2 [1 −
621
to investigate this, this process was excited with a sinusoid of amplitude 0.1 m and 0.05 cycle/sample or 20 samples/cycle signal for various operating point. The results are shown in Figure 7. It is clear from the figure that nonlinearity increases as the operating point moves from the center of the tank, i.e., 0.5 m. This is also expected because the curvature of the inner surface of the tank increases as it moves away from the center.
is chosen at 0.4m instead of 0.5m because the process is less nonlinear in the vicinity of the latter operating point. 5.1 Dependency of nonlinearity on the amplitude and frequency of input signal All linear system exhibits ‘sinusoidal fidelity’. Sinusoidal fidelity states that a sinusoidal input to a linear system will produce a sinusoidal output of identical frequency. The amplitude and phase of the output sinusoid may differ from the input sinusoid. But this is not the case for a nonlinear system. Generally, in response to a sinusoidal excitation a nonlinear system will produce new frequencies in addition to the original frequency of the input sinusoid. This property requires the use of a sinusoid as an excitation signal for testing nonlinearity of a system. In this example, a sinusoidal excitation (added as a deviation to the steadystate Fi ) with varying frequencies and magnitudes are used. The amplitudes vary from 0.01 to 0.2 and the frequencies vary from 0.015 to 0.1 cycles/sample or 54 smaples/cycle to 10 samples/cycle. 4096 data points were used for each nonlinear analysis. The results are shown in Figure 6. This figure shows that for small magnitude input signals, the process can be assumed locally linear. As the amplitude of excitation goes up, so does the nonlinearity. Change of frequency of the input signal is not very sensitive to the nonlinearity analysis. Therefore, it can be concluded that nonlinearity of the process strongly depends on the size of the excitation or input signal.
6. CONCLUSIONS A method for better estimation of bicoherence has been presented. The modified N GI and N LI perform better in detecting signal non-Gaussianity and nonlinearity. A new index, Total Nonlinearity Index (T N LI), has been developed and its performance in quantifying signal nonlinearity is demonstrated using illustrative simulation examples. These indices are applicable to any time series signal and does not require a sinusoidal excitation of a system. Further research is in progress towards the use of these bicoherence based indices for root cause diagnosis of plant-wide oscillations. REFERENCES
TNLI
4 3 2 1 0 0.1
0.2 0.08
0.06 0.04 w, frequency 0.02
0.15 0.05
0.1 A, amplitude
0 0
Fig. 6. Dependency of nonlinearity on the input signal’s amplitude and frequency 3.5 3
TNLI
2.5 2 1.5 1 0.5 0
0.2
0.4 0.6 Nominal Level
0.8
1
Fig. 7. Dependency of nonlinearity on the operating region of the system 5.2 Operating point Dependency of nonlinearity It is well-known that the nonlinearity of a system depends on the operating point of the system. In order
622
Agrawal, P. and S. Lakshminarayanan (2003). Tuning pid controllers using achievable performance indices. Ind. Eng. Chem. Res. 42, 5576–5582. Choudhury, M. A. A. S. (2004). Detection and diagnosis of control loop nonlinearities, valve stiction and data compression. PhD thesis. University of Alberta, Canada. Choudhury, M. A. A. S., S. L. Shah and N. F. Thornhill (2004). Diagnosis of poor control loop performance using higher order statistics. Automatica 40(10), 1719–1728. Collis, W.B., P.R. White and J.K. Hammond (1998). Higher-order spectra: The bispectrum and trispectrum. Mechanical Systems and Signal Processing 12, 375–394. Fackrell, J.W.A., P.R. White, J.K. Hammond, R.J. Pinnigton and A. T. Parsons (1995). The interpretation of the bispectra of vibration signals - i. theory. Mechanical Systems and Signal Processing 9, 257–266. Helbig, A., W. Marquardt and F. Allg¨ower (2000). Nonlinearity measure: definition, computation and applications. Journal of Process Control 10, 113123. Hinich, M. J. (1982). Testing for gaussianity and linearity of a stationary time series. Journal of time series analysis 3, 169–176. Nikias, C. L. and A. P. Petropulu (1993). HigherOrder Spectra: A nonlinear signal processing frame work. Prentice Hall. New Jeresy. Rao, T. S. and M. M. Gabr (1980). A test for linearity and stationarity of time series. Journal of time series analysis 1(1), 145–158. Thornhill, N. F., F. Cox and M. Paulonis (2003). Diagnosis of plant-wide oscillation through data-driven analysis and process understanding. Control Engng. Prac. 11(12), 1481–1490.