Comparison of Nonlinearity Measures Based on Time Series Analysis for Nonlinearity Detection

The main purpose of this paper is a study of the efficiency of different nonlinearity detection methods based on time-series data from a dynamic process as a part of system identification. A very useful concept in measuring the nonlinearity is the definition of a suitable index to measure any deviation from linearity. To analyze the properties of such an index, the observed time series is assumed to be the output of Volterra series driven by a Gaussian input. After reviewing these methods, some modifications and new indices are proposed, and a benchmark simulation study is made. Correlation analysis, harmonic analysis and higher order spectrum analysis are selected methods to be investigated in our simulations. Each method has been validated with its own advantages and disadvantages.


Introduction
System identification involves obtaining the structure and the unknown parameters of the model of a system from the input-output information.If the structure of the equations is unknown, the problem is called structure identification.For the identification, the input/output information is required to be available from suitable tests conducted on the system.
The input and output might have linear or nonlinear dynamic relation.Linear identification is much simpler than nonlinear identification, so it is the first choice in system identification.However, in many applications, the effect of nonlinearity cannot be ignored.In such cases nonlinear system identification has to be carried out Choudhury et al. (2008).
Before proceeding with any testing and identification, it is useful to establish whether the system can be considered linear across the desired operating range.It is important to be able to measure the degree of non-linearity of a process under various input excitation signals or operating conditions.
In the literature, many tests have been proposed for measuring the nonlinearity but they do not provide a general answer to the problem, since they rely on different assumptions.A simple method is injecting a series of single sine waves of increasing amplitudes and look for distortion of the output signal Billings and Fadzil (1985).
Testing linearity of time series can be done using parametric or nonparametric methods.In a parametric approach, the nonlinear distortion depends on an unknown vector parameter θ, and it will be zero when the vector parameter (θ) is zero Saikkonen and Luukkonen (1988).
The advantage of the parametric approach is that relatively few observations are required.The disadvantage is that the tests are generally sensitive to a misspecified nonlinear model.The nonparametric methods are based on analysis of input/output data avail-able, or just output measurements.
There are some papers for nonparametric methods in the Time-domain Hjellvik and Tjøstheim (1995); Barnett and Wolff (2005) but most of these methods have been applied in the frequency domain.In fact, the nonparametric approach has been based mainly on higher order spectra.Initially, Rao and Gabr (1980) implemented Brillinger (1965) method for measuring the deviation of a process from linearity and Gaussianity by estimation of the bispectrum of observed data.Then Hinich (1982); Ashley et al. (1986), modified it by introducing a chi-squared statistic index.
The advantage of the non-parametric test compared to parametric methods, is the invariance with respect to linear filtering of the original data.A disadvantage is that more observations are necessary to obtain a power comparable to that of the best parametric test for a fully specified alternative Rusticelli et al. (2008).Caillec and Garello (2004) have made a comprehensive comparison between the above methods to validate their efficiency.
All of these papers are based on bispectrum analysis, while there are some other nonparametric methods which have been presented in the literature of nonlinearity test: Linear spectral density method in Bendat and Piersol (1980), linear correlation method in Szücs et al. (1975), nonlinear cross correlation methods in Billings and Voon (1986), and higher order auto correlation method in Billings and Voon (1983).
In order to investigate these methods, Haber (1985) has made a survey of the different methods and compares them.There are also reported researches on test for bilinearity in Varlaki et al. (1985) which are out of the scope of this paper.The main novel contributions of the present paper are as follows: 1) A comprehensive mathematical description of the nonparametric approach with the emphasis on frequency domain methods is done.In the different sections, we derive properties of these indices under various conditions.The efficiency of the different methods for detecting the nonlinearities is reported.An investigation has been performed with respect to additive white measurement noise.
2) By the mathematical description of each method and investigating their indices, some new indices are presented.
3) With the focus on higher order spectra and harmonic analysis methods, the higher order cross correlation function has been reviewed.Also some new indices based on the harmonic analysis method have been proposed.
4) A benchmark simulation study is made for selected methods.A comparison among them has been done in order to investigate their efficiency and appli-cation especially from a practical viewpoint.
Single input-single output (SISO) systems, time invariance, and asymptotically stability of the process are three basic assumptions which are made throughout this paper.The outline of the paper is as follows.In section 2, we review and compare several nonparametric methods associated with nonlinearity detection indices.We survey mathematically the higher order auto-and cross-spectra based methods with the emphasis on the type of the input signals in more details in section 3. Some simulation results in order to test these methods and to investigate their efficiency regarding to measurement noise and different types of inputs are presented in section 4.
Bicoherence of y

Nonlinearity measure methods
This section discusses various measures of the degree of nonlinearity.We start our study by investigation of some nonlinearity measures which are presented in Haber (1985).All methods for measuring the nonlinearity, which are presented here, are based on the block diagram in Figure 1, i.e. open loop nonlinear process with additive measurement noise.
The basic principle behind these methods is using some analytical properties of LTI (linear time invariant) systems.Then the nonlinearity is measured based on deviation from linearity.The inputs which are  (1918): where the notation ΣΠ means summing over all distinct ways of partitioning η 1 η 2 . . .η 2N into pairs.Signals such as sine wave, pseudo random binary signal (PRBS) conform to this class of inputs Billings and Voon (1986).However if a Gaussian input passes through a nonlinear static block, this theorem can still be used as a nonlinearity test, but for non-Gaussian input there is no analogous theorem Nicholas et al. (2009).

Linear cross correlation (frequency domain)
For LTI systems, the principle of superposition holds and the system is appropriately described by one or more independent modes that are invariant to the level of excitation.So, the auto-and cross-spectral densities associated with the system input/output data can be shown to be sufficient to characterize (identify) the system.This analysis, however, only considers second-order moments of the input/output data Ljung (1999).It has been shown for any LTI systems without noise, the squared coherence function between input and output is one.(see Bendat and Piersol (1980)): Applying the coherence function based nonlinearity indices with a random variable as input, may give us correct decision, but it does not have this property with periodic deterministic signals as the excitation.If u(t) is any periodic function of time, measured over N periods, and there is no noise on u(t) and y(t), then the values of φ yu (jω), φ uu (jω) and φ yy (jω) are independent of N.So the effect of nonlinearity will be ignored when a periodic signal is used McCormack et al. (1994).Therefore, if the LTI system is driven by a random excitation input and there is no noise then the following index is zero Haber (1985): However this index is affected by additive noise.Because of noise, unknown disturbances, and biases and variances in estimation of the coherence function from a finite data set, the maximum value of coherence function may give us incorrect decision.Since NLI1 selects the minimum value, it causes high sensitivity results.Hence the NLI2 index is proposed.
By the definition of nonlinearity (deviation from linearity), this index captures the deviation at worst case condition.Meanwhile, the average of squared coherence function may obtain more robust performance.This motivates the definition of the N LI3 index.
Lemma 1: If the system is LTI, then: If the system is nonlinear or there is additive white noise, then N LI1, N LI2, N LI3 > 0.
Proof: See Appendix.
In section 4, these indices are evaluated and compared to each other.

Linear cross correlation (Time-domain)
Based on correlation analysis methods in the frequency domain, two discrete and continuous nonlinearity indices can be defined in the Time-domain.With the assumption of no output noise and stationary input signal, it is proved in Haber (1985) that if the system is linear, then: is zero.Analogously, this index is also affected by additive white noise.

Nonlinear cross correlation method
Assume a Gaussian white noise or pseudo random input signal is applied to the system, and define a normalized form of the squared input as: . Then it is proved that the cross correlation function φ yx (τ ) is zero for linear systems and differs from zero for nonlinear systems Haber (1979): This leads to a nonlinearity index in the discrete Time-domain: The nonlinearity index N LI5 > 0 is robust to additive white noise as its effect will vanish asymptotically Haber (1985).

Higher order auto correlation method
By using a white Gaussian random or pseudo random input signal of zero mean for measuring the nonlinearity, it is shown in Billings and Voon (1983) that the higher order autocorrelation function φ yy 2 (τ ) of the normalized output signal is zero for linear systems.Let In order to get the correlation values between -1and 1, we normalize the square of the output signal as well This leads to the following index in the discrete Time-domain: Analogously, the nonlinearity index N LI6 > 0 is robust to additive white noise Haber (1985).If any system just includes odd order (polynomial) nonlinearity, the indices will be zero.So they can detect odd-order nonlinearities only if the mean value of the input signal is not zero but the even-order nonlinearities can be detected without such condition.

Harmonic analysis
Assume u (t) = Asinωt is applied to a system.Since each periodic signal has a Fourier series, the output signal can be represented as Since any LTI system generates no more harmonics than found in the input signal, this leads to a nonlinearity index as below Haber (1985): where If the system is LTI and there is no noise, then N LI7 = 0.

Harmonic matrix
As mentioned earlier, when the system is nonlinear, some extra harmonics will be generated in the output signal.Based on this property, some new nonlinearity indices are proposed.Assume there is a Fourier series representation of u(t) and y(t).So: Consider the following integrals: Here N, M is the number of harmonics of input/output signals, respectively.So based on Cauchy-Schwarz inequality: where c i and d j are Fourier coefficients (harmonics) of input and output signals respectively.In this index, it is true that M ≥ N .If the system is linear then M = N → N LI = 1 − Γ = 0, when the system is nonlinear M > N , so some extra terms appear in denominator and N LI > 0. Based on above discussion, it is possible to construct the relation between input and output harmonics as a matrix: By solving the above equation and finding the H matrix we can deduce that: If H is diagonal then the system is linear.If H is non-diagonal then the system is nonlinear.So by defining the following index we can detect the nonlinearity of the system.

Harmonic distortion
Consider the static nonlinear system y = 3 p=1 u p excited with a multi-harmonic periodic signal.The output signal includes linear, quadratic and cubic nonlinearity.The frequency spectrum of the output has the property that a nonlinear system generates more harmonics than the basic harmonics.
Quadratic nonlinearity creates additional harmonics without main harmonics, while cubic nonlinearity generates both basic and new harmonics in the output Pintelon and Schoukens (2001).Here a multi-harmonic signal is defined as: with double sided spectrum: where A, i, φ are vector of amplitudes, harmonic and phase, respectively, and F is the number of harmonics.
It is noticed that if the phase of the input signals is close to zero, the amplitude of new harmonics become maximum but there may be a problem due to maximization of excited signal in Time-domain if it causes deviation from the current operating point.So, if we can manipulate the input signal, the phase should be selected so that the CF (crest factor) of it gets minimum Solomou and Rees (2003).
In general form, consider a family of static nonlinear systems as a p u p The output response for a multi-harmonic signal is shown in Evans et al. (1994) as: It is easily found that if there are 2F harmonics in input signal, then at most (2F ) p harmonics will be produced in the output signal, where p is the nonlinearity order.Assume only a cubic nonlinearity, so the output will consist of 8 terms with three terms at ω 0 , one term at 3ω 0 , three terms at −ω 0 and one term at −3ω 0 .Analogously with just a quadratic term there would be 4 terms which two terms are at 0, one term at 2ω 0 , one term at −2ω 0 .Usually the terms (contributions) generated by the odd nonlinearity are considered as two types.
Type I Contributions: These contributions are generated by combinations of a test frequency with a pair of equal positive and negative frequencies, which produces a contribution at the test frequency itself.For example, the combinations δ (ω − (ω 0 − ω 0 + ω 0 )) and δ (ω − (ω 0 + 3ω 0 − 3ω 0 )) will each produce a contribution at ω 0 .
Type II Contributions: They are generated by combinations of signal content frequencies which are not included test frequency but result in test frequency.
For example, the combinations δ (ω − (3ω 0 − 9ω 0 + 7ω 0 )) and δ (ω − (9ω 0 − 3ω 0 − 5ω 0 )) will each produce a contribution at ω 0 .It is possible to design signals which do not suffer from any type II distortions by proper selection of their component harmonics.This is called NID (no inter-harmonic distortion) Evans et al. (1994).Assume systems of the form eq.( 24): Lemma 2: For a F-harmonic excitation signal with a same amplitude A and excited frequency ω 0 , the number of harmonics (N p ω0 ) and the number of new harmonics (N p new ) are: Proof: See Appendix.Lemma 3: For the NID input excitation signal, the number of Type I harmonics are: Referring to Lemmas 2-3, the number of total harmonics of a cubic nonlinearity in test frequency (ω 0 ) are summation of Type I and Type II contributions . Now a nonlinearity index is required to measure the power of produced additional output components (harmonics).The squared coherence function could be a suitable indicator Evans et al. (1995a).Hence an index is defined as where Ω k is set of total frequencies in input signal and ω e is an excluded frequency (new harmonics set).For a stable LTI system, since no additional harmonic is generated, we have Y (jω e ) = 0 → N LI9 = 0 otherwise N LI9 > 0. Another suitable index which can be defined is This index measures the ratio of power of new harmonics to total harmonics.These indices are influenced when considering additive white noise.Beside the characteristics of additional (new) harmonics, the number of new harmonics is also important and it seems well if they contribute to the nonlinearity index.Hence a new index which includes both the number and amplitude of harmonics for measuring the nonlinearity is presented as below.
where α is a parameter which can be selected by practitioner.It should be noticed that additional harmonics may appear at frequencies higher than the system bandwidth.However, if the system under test has much smaller bandwidth than the wide band excitation signal, then these high frequency harmonics may be attenuated.

Higher order nonlinearity measures
This section focuses on nonlinearity indices based on the higher order moments.The third order moment can be viewed as a powerful tool for measuring the nonlinearity, particularly in time series analysis.This method has been used in several references Marzocca et al. (2008); Nicholas et al. (2009).As it could be seen here, the bispectrum analysis has some advantages.First, its value for LTI systems is independent of bifrequencies.Second, many statistical tests such as Hinich (1982) can be defined for it.Third, it is able to capture quadratic nonlinearity with high efficiency.Fourth, its estimator converges asymptotically, although, it needs a large amount of data.In the following, we review and study the properties of the Volterra series as a mathematical representation for nonlinear systems study.
The Volterra series is a mathematical tool widely employed for representing the input/output relationship of nonlinear dynamic systems.Volterra series are based on the expansion of the nonlinear operator representing the system into a series of integral operators.These integral relationships are completely defined by the Volterra kernels.
An alternative representation, which generates the frequency domain counterparts of the Volterra kernels, can be viewed as a generalization of the usual frequency response function for nonlinear systems.A fairly large class of dynamic systems can be modeled according to these concepts and are accordingly represented in terms of Volterra series Nijmeijer and Schaft (1990).The system output response can be obtained as a sum of individual components: where Such expressions state that the output y(t) is related to the input u(t) through multi-dimensional convolution integrals of the input with the functions h i .It is worth remarking that h n (τ 1 , τ 2 , . . ., τ n ) = 0 ∀τ i < 0, in order for the system to be causal.Existence and uniqueness of the kernels, as well as convergence analysis of the series can be found with relevant details in Schetzen (1980); Evans et al. (1995b).Volterra series in the frequency domain can be written as .e −iω1τ1 e −iω2τ2 . . .e −iωnτn dτ 1 dτ 2 . . .dτ n (37) As in the LTI case, these complex functions carry the same information about the system as their Timedomain counterparts, but they are often easier to work with, both experimentally and analytically.The output property of a system represented by Volterra series in the presence of different types of inputs is presented in Bedrosian and Rice (1971).

Time-domain nonlinearity test
Consider Figure 1 and assume that all the odd-order moments of the input signal u(t) are zero and the noise e(t) is an independent, zero mean stationary random signal.It has been shown that if a system is LTI, then et al. (2009).Hence, whenever φ zz 2 (τ ) = 0, this indicates that the system under test may be nonlinear.
Note that the test distinguishes between linear additive noise corruption of the measurements and distortion due to nonlinear effects.Based on this function a nonlinearity index can be defined as The advantages of this index are its use of higher order correlations and its robustness to additive white noise.However, the disadvantages could be named as being limited to Gaussian inputs and lack of normalization.
There is another form of higher correlation Timedomain method which can be used for measuring the nonlinearity.It is based on cross correlation of the input and output signals.Keeping the previous assumptions, the process under test is LTI if φ zuu (σ 1 , σ 2 ) = 0 ∀σ 1 , σ 2 (Billings and Chen (1980); Billings and Voon (1986)).This is also true if the input belongs to the separable class of random processes Nuttal (1958); Collis et al. (1998).Based on this function another index can be defined as It has the same advantages and disadvantages as N LI12.

Frequency domain nonlinearity tests
The analytical expressions for the higher order spectra are investigated here using a Volterra series approach on the assumption of a zero-mean, stationary Gaussian input.The bispectrum is the lowest order of the higher order statistics (HOS) and has been shown to be a promising tool for both detecting and identifying system nonlinearities.
The statistical properties associated with bispectrum estimation and theoretical investigations involving bispectrum estimates were reported in several references (for instance Caillec and Garello (2004)).
Based on the formula for the first and second order Volterra kernels, an analytical expression for the output bispectrum φ yyy can be derived.However, a closed-form solution for the auto-bispectrum for a general class of nonlinear systems has not yet been presented.
It has been shown that φ yyy can be expressed as a function of the first and second order Volterra kernels H 1 and H 2 and the input spectrum φ uu Marzocca et al. (2008).Higher order auto and cross spectra can be defined based on the third order moments or cumulant.The moments and cumulants of a stochastic signal y(t) are defined as Fackrell (1997): Bispectrum is the frequency domain counterpart of the third order cumulant and is defined for a stationary signal y(t) Since for a zero mean stationary signal c 3 y = m 3 y , the bispectrum can be written as Priestley ( 1988) Thus, an analytical expression for the bispectrum is obtained by substituting the Volterra model into the above equation.In a simple form of nonlinear system which is modeled by two Volterra kernels (H 1 , H 2 ) and using eq.( 36) and eq.( 43), a closed formula for bispectrum has been presented in Marzocca et al. (2008).The bispectrum contains the H 2 term, only making this quantity suitable to capture quadratic nonlinearities.When cubic or higher order nonlinearities are to be detected, HOS, e.g. the trispectrum should be utilized.
As nonlinear terms appear in the classical linear spectral analysis as a distortion, in an analogous way, nonlinear components of higher order than the quadratic ones may appear also in the bispectrum.Analogously, the higher order cross correlation between u(t) and y(t) are defined as: By the generalization of all Volterra kernels and substituting in eq.( 45), the odd-order kernels are eliminated and even-order kernels may be expressed as integral functions of H 4 , H 6 , and the power spectrum of input.It has been shown for Gaussian white noise input with constant auto spectrum (P) and for n evenorder kernels Worden and Tomlinson (2001) dω (1) . . .dω ( n 2 −1) (46) The most important consequence of eq.( 46) is So as P → 0, the higher order cross correlation function tends to H 2 .

Higher order spectra as nonlinear index
Regarding to theorem (1) and eq.( 43), the bispectrum will produce nonzero values, when the input to a linear system is non-Gaussian.So the appearance of nonzero values is not only due to the nonlinearity, but may also be due to the non-Gaussian nature of the input.So for a LTI system defined by its transfer function (H 1 ) which is subjected to a zero-mean Gaussian input, an expression for the output bispectrum is readily obtained as Since Bis u (ω 1 , ω 2 ) = 0, then Bis y (ω 1 , ω 2 ) = 0.This leads to the nonlinearity index: Bis y (ω 1 , ω 2 ) (49) Its advantages are using higher order spectra and robust to additive white noise while weaknesses are limited to Gaussian input signals and not normalization.Based on eq.( 45), it is easily shown that for a Gaussian input signal and LTI system φ yuu (ω 1 , ω 2 ) = 0∀ω 1 , ω 2 which can be used as a nonlinearity index.This function has the same advantages and disadvantages as the bispectrum.A convenient normalization of the bispectrum is called bicoherence and is obtained by dividing the magnitude of Bis y by the appropriate product of output spectra Priestley (1988).
where P y (ω 1 ) is power spectrum of signal.Using the Cauchy-Schwarz inequality, the normalized version of bispectrum is as below: For a linear system, the bicoherence function is constant as where If the input is Gaussian then µ 3 = 0 and bic y = 0.The property of bicoherence constancy of an LTI system to any bifrequencies, allows some indices to be defined for white Gaussian input.Flatness test is an index which has presented in Choudhury et al. (2008): The advantages of this index are its use of higher order correlations, and normalization while its disadvantage is limitation to Gaussian inputs.The N LI15 has another disabililty when the bicoherence function generates similar and large non zero quantities.In this condition, using N LI15 gives low value so it is concluded the system is linear while it is nonlinear.So we define another nonlinearity index as: With the assumption of non-Gaussian signal, the following indices can also be defined to detect the nonlinearity.
Similar to linear cross correlation method, these indices have their own advantages and disadvantages.Analogously, a convenient normalization of φ yuu is obtained as Using the Cauchy-Schwarz inequality, the normalized version of φ N yuu satisfies: Based on the normalized higher order cross correlation function, the following indices can be defined These indices have the same properties of the previous ones.For all above indices, nonlinearity can be tested by N LI > T where T is a given threshold that should be determined by the practitioner.

Higher order spectra with additive noise
In this section we review and investigate the effect on higher order spectra caused by additive white Gaussian noise on the output ( z = y +n).Due to an assumption of Gaussian white additive noise, it could be shown that the higher order auto and cross spectrum functions are robust to white noise.However the normalized versions of these functions, because of appearance of some additional terms in denominator, are affected with white Gaussian additive output noise.Using the following lemma causes the increase of the accuracy of nonlinearity detection with noise.Lemma 4: For a system with Gaussian input and additive white noise, the normalized higher order auto and cross correlation functions satisfy eq.( 84) and eq.( 90).Proof: See Appendix So, the following indices can be defined: By considering this lemma we can decide more accurately about the nonlinearity of the system provided that we know the signal to noise ratio.

Simulation results
In this section, we compare the advantages and disadvantages of the most promising nonlinearity indices using a benchmark simulation example (Duffing system).
The Duffing system is a well-known nonlinear second-order differential equation used to describe many physical and engineering problems.The equation is given by m d 2 dt 2 y (t) + c d dt y (t) + k 1 y (t) + k 2 y 2 (t) + k 3 y 3 (t) = u(t) where m, c, k 1 , k 2 , k 3 and u are mass, damping coefficient, linear stiffness, nonlinear (quadratic) stiffness, nonlinear (cubic) stiffness and excitation input, respectively.
In the conventional Duffing equation, the quadratic term is set to zero and it is called symmetric Duffing oscillator.However, many successful applications in engineering have been reported by expanding it to asymmetric one (k 2 = 0).This system is time-invariant and can be unstable at high levels of excitation in the cases when k 3 ≤ 0.

Nonlinearity measure based on coherence function
In this section, we evaluate the performance of coherence method based nonlinearity indices.Assume a periodic function as u (t) = Asin(ω 0 t) is the excitation signal of the Duffing system without any noise.The results in Figure 2 are obtained by recalling the coherence function formula in eq.( 1) and three indices N LI1, N LI2, N LI3 together with the estimation of γ 2 yu .
In this case the coherence function is estimated for four parameter settings of the Duffing system as linear (k 2 = k 3 = 0), low order, medium and high order nonlinear (according to different values of k 2 , k 3 ).For each of the four parameter settings, the nonlinearity degree based on three indices N LI1, N LI2, N LI3 is zero.It is confirmed that the coherence-based nonlinearity indices are zero when a periodic excitation signal is used, and the indices fail.By simulating the previous scenarios, but using non-periodic excitation signal (PRBS) we have the results in Figure 3.
Figure 4 shows the same experiments but with additive output noise.The index N LI1 is not sensitive to nonlinearity and always gives a very close value to zero.Although N LI2 gives the maximum defection of nonlinearity, it is overly sensitive.The N LI3 detects more correctly the level of the various nonlinearities and noise with different variance.In conclusion, the following notes are made concerning the coherence function method: 1-The coherence function is unable to detect the nonlinearity with periodic input signals.2-It is affected by additive white noise (setting a threshold depends on noise variance).3-The N LI3 is better than N LI1, N LI2.It has a tradeoff between capturing the nonlinearity and sensitivity to noise, since N LI3 = 0 for the linear system due to noise.4-Because of the finite data in spectrum estimation, the N LI1 and N LI2 fail even in noise free cases.

Nonlinearity measure based on harmonic analysis
In this case, the efficiency of the harmonic analysis approaches has been evaluated.Using multi-harmonic input u(t) as an excitation signal, we get the results shown in Figure 5 in the case of no output noise.By using the multi harmonic signal, more extra harmonics are generated, so it is more able to distinguish between linearity and nonlinearity by both indices.
It should be noticed that because of the numerical simulation of nonlinear system, some noise behaviors in the vicinity of each given frequency will be present in the harmonic plot.
Figure 5 illustrates that the measurement noise affects some of the new harmonics in the nonlinear system, but it is still possible to catch the nonlinearity by other new harmonics.However this figure shows  that the two nonlinearity indices in N LI9, N LI10 suffer from additive Gaussian noise, and may end to wrong decision.In summary, the following conclusions are made: 1-The harmonic analysis method is affected by noise output.
2-Not all harmonics in the excitation signal are effective in this method.
In fact those harmonics which are designed based on the bandwidth of the system are functional.Although noise deteriorates the ability of this method, because of its high frequency feature, the method is more robust if the applied harmonics are inside the system bandwidth and have values in low frequencies.

Higher order spectrum analysis
In this section, the efficiency and functionality of the higher order spectrum methods have been investigated.Assume the Duffing system is driven by a PRBS sig-nal.The amplitude of the excitation signal should be selected correctly.If it is too low the nonlinearity is not detectable and if it is too high it can sometimes lead to wrong decisions.
In this simulation the performance of four nonlinearity indices N LI15, N LI16, N LI19, N LI20 considering additive white noise meantime are compared.Figure 6 shows that for a LTI system with Gaussian input, the two normalized higher order auto-and cross-correlation functions are zero and for nonlinear mode are non-zero.The effect of additive Gaussian noise has been shown in Figure 7.By comparing Figures 6-7, the indices NLI16, NLI20 which find the maximum value of higher order spectrum functions work better because of the unknown behavior of the noise spectrum that may cause the bifrequencies functions tend to more flatness.It can be seen in the figures that this index is better when noise exists.Here, by reconstructing the higher order auto and cross correlation functions through eq.( 84) and eq.( 90), in appendix, the following results in Figure 8 are obtained.It is illustrated, in Figure 8, that the nonlinearity indices N LI23, N LI24 are modified from the N LI16, N LI20 indices while the measurement noise exists.In conclusion, the following notes are made: 1-Both the auto (bicoherence) and higher order cross spectrum functions could be effectively used for nonlinear detection.
2-Bicoherence is used for time series analysis where just output measurement is available.
So for problems such as modeling and identification, the cross correlation method seems better to use.The indices N LI16, N LI20 are better than N LI15, N LI19 indices.

Conclusion
In this paper, we have reviewed and compared the efficiency and performance of several time series based nonlinearity measures for nonlinearity detection.This study is based on the comparison of their functionality, sensitivity to noise, complexity, dependency on type of input and efficiency in capturing the nonlinearities.
Four methods are investigated and for each one some new nonlinearity indices have been introduced.Since the higher order statistics analyzes are stronger, the higher order spectra (auto and cross) are studied in more details.Relevant to each method some useful characteristics and modifications are presented.
A nonlinear benchmark system which is used in many researches (Duffing system) has been used in selected methods to study their efficiency.Particular attention is given to the effect of output noise.
In conclusion, a summary of several nonlinearity indices properties has been presented as table 1.In this table the symbol (*) remark new nonlinearity indices and symbols (-, +, ++) show mediocre, good and better efficiency respectively..δ(ω− For simplicity and with no loss of generality, we can assume that the amplitude of multi-harmonic signal is fixed to A and phase is zero.So: It is true that the number of harmonics relevant to ω 0 in the output signal is proportional to: Regarding to eq.( 72), the total number of permutations of different values of i (n 1 ) , i (n 2 ) , i (n 3 ) so that i m = i (n 1 ) + i (n 2 ) + i (n 3 ) = 1, ∀m = 1 . . .2F, and for F = 1, 2, 3, 4 are shown in table 2. In fact this rule is done in 2F steps by fixing i (n 1 ) to its values and find the permutations of i (n 2 ) , i (n 3 ) so that i m = 1.
Following the above rule in table 3, It is deduced that N 2 ω0 = 2F − 2. Moreover the number of zero harmonics in the output signal is needed.So using previous procedure generates table 4. By figuring out the result of table 4, the number of zero harmonic can give as N 0 = 2 Clearly, for F harmonic in excitation signal, there are (2F ) 3 total harmonics.So the number of new non zero harmonics In quadratic nonlinearity, it is readily finding that N 0 = 2F .Therefore, the number of new non zero harmonics are and the number of positive new harmonics is

Figure 2 :
Figure 2: Coherence function of Duffing system in linear and nonlinear modes with a periodic excitation signal

Figure 3 :Figure 4 :
Figure 3: Top.Coherence function of Duffing system in linear and nonlinear modes with nonperiodic excitation signal.Bottom.Nonlinearity measures NLI1, NLI2, NLI3 for Duffing system with non-periodic excitation signal

Figure 5 :
Figure 5: New harmonics generation of Duffing system in linear and nonlinear modes with a multi harmonic excitation signal and Top.without measurement noise.Bottom.with measurement noise

Figure 6 :
Figure 6: Nonlinearity measure based on indices N LI15, N LI16, N LI19, N LI20 for Duffing system in Top.linear mode.Bottom.nonlinear mode

Figure 7 :Figure 8 :
Figure 7: Nonlinearity measure based on indices N LI15, N LI16, N LI19, N LI20 for Duffing system in nonlinear mode with measurement noise

2 .
Proof of Lemma 3:Since the NID signal does not include any Type II contributions, on the other hand, Type I contributions are generated by combinations of a test frequency with other pair of equal positive and negative

Table 1 :
Summary of different nonlinearity indices

Table 2 :
Number of permutations of harmonic ω 0 (N 3 i 1 i 2 i 3 i 4 i 5 i 6 i 7 i 8

Table 4 :
Number of permutations of zero harmonic (N 0 ) i 1 i 2 i 3 i 4 i 5 i 6 i 7 i 8

Table 5 :
Number of Type I harmonics (N 3 ω0