A Digital Method of Evaluating Nonlinear Frequency Response Functions Based on Conditioned Spectral Analyses

This paper focuses on the issue of evaluating nonlinear frequency response functions (NOFRFs) of the nonlinear systems with a general input. A new digital approach of evaluating NOFRFs is proposed based on the method of conditioned spectral analysis (CSA). Firstly, a multiple-input/single-output (MISO) linear system with a series of power characterized inputs is obtained by decomposing a nonlinear system according to the NOFRFs theory. Secondly, the correlations among the inputs of various orders are removed by applying CSA method, obtaining an algorithm of identifying the NOFRFs and evaluating the contributions of different order nonlinearities to the output of the system. In the CSA procedure, the different order inputs are conditioned in the sequence from the first order to the highest order. Lastly, two kinds of nonlinear systems with a general input are simulated to verify the effectiveness of the proposed approach on evaluating the frequency response properties of nonlinear systems. It is shown that the results reached by the proposed method are very close to the numerical results obtained by the fourth order Runge-Kutta method, verifying that the proposed method is very effective on the evaluation of NOFRFs.


Introduction
Most practical mechanical systems have nonlinear dynamic characteristics. When a general input is acted on an unknown nonlinear system, the system's nonlinear dynamic behaviors can generate multiplication and intermodulation frequency components in the output spectra. Due to the fact that the higher order frequencies are much richer than the input frequencies and multiple interferences may occur among all of the frequency components, difficulties are often encountered when one wants to evaluate the nonlinear system's frequency response functions (FRFs) or to quantify the effects of the input on the system. A common method of analyzing the frequency response properties of nonlinear systems is based on Volterra series [1], which provides a general framework by introducing generalized frequency response functions (GFRFs) [2,3]. Because the GFRFs represent the nonlinear frequency response features of various orders, the methods of identifying GFRFs has been researched in many studies. One of the identification methods was implemented by minimizing the mean square deviation of output under the condition that the input is a zero-mean stationary Gaussian process [4][5]. But in many mechanical systems, the input does not satisfy Gaussian statistics, such as the vibration excitations generated by raceway geometrical errors in rolling bearings [6]. Therefore, an extended digital method was proposed based on higher order spectra for estimating the linear and quadratic transfer functions of the quadratic time-invariant system with one general input [7,8] and two general inputs [9]. However, for the reason that the number of frequency interactions increases quickly with the order of GFRFs [10][11][12], explicit expressions of frequency interactions required in the extended approach become very complicated when the order is higher than two or three. Another way to evaluate the GFRFs was by applying parametric modeling techniques, typically involving orthonormal basis expansion [13], analytical derivation [14], and harmonic probing [15][16][17][18][19] which has been applied to detecting cracks in rotor-bearing systems and cantilever beams [20]. Although the methods of para-metric modeling can obtain any order GFRF theoretically, they require not only the output data available but also the knowledge of governing motion equations of system, and inherently have the problem of convergence in iteration processes [21,22]. Recently, a new concept of nonlinear output frequency response functions (NOFRFs) has been proposed by Lang and Billings [23], which allows the analyses of nonlinear systems to be carried out in a manner similar to linear systems and provides a straight forward explanation for the output frequency response properties of nonlinear systems. The effectiveness of NOFRFs has been extensively verified by Peng and Lang etc., involving energy transfer of bilinear oscillators [24], comparison between NOFRFs and harmonic balance of Duffing oscillators [25], parameter estimation of multi-degree-of-freedom (MDOF) nonlinear systems [26] and applications in structural damage detection of nonlinear systems [27] etc. In these studies, however, the NOFRFs are obtained by the algorithm of least square, which requires the system be excited N times with different input signals and the number of excitations be equal or greater than the highest order of NOFRFs, and thus the algorithm is subject to the constraint that the system input should be controllable either in simulations or in practical experiments. In this paper, an alternative approach is introduced to evaluate the NOFRFs of the stable time-invariant nonlinear system with a general input. This method is performed by using the technique of conditioned spectral analysis (CSA) which can remove all of the interferences among different order output frequency responses. An algorithm is provided and the contributions of each order nonlinearity to the output are distinguished without requiring that the system be excited by different input signals many times.

NOFRFs of Nonlinear Systems
For a stable time-invariant nonlinear system with a general input and a noise contaminated output, the input-output relationship in the neighborhood of an equilibrium point can be expressed by Volterra series as (1) in which (•) and (•) represent the system input and output respectively, ( ) denotes the output components generated by the nth order nonlinearity ( = 1, ⋯ , ), ℎ ( 1 , ⋯ , ) is called 'the nth order Volterra kernel' or 'the nth order impulse response function' which can reduce to the well-known impulse response function of linear systems when = 1, ( ) denotes the noise imposed on the output, represents the maximum order of nonlinearities considered for the system. The output components of the orders higher than are omitted. In frequency domain, the input-output maps of Eq. (1) can be formulated as [11]   where signifies the angular frequency, i.e., = 2 , (j ) and (j ) represent the system input and output spectra respectively, (j ) denotes the nth order output spectral components, d refers to the area of a minute element on the n-dimensional hyperplane = 1 + ⋯ + , (j 1 , … , j ) represents the Generalized Frequency Response functions (GFRFs) equaling to the multidimensional Fourier transform of the nth order impulse response ℎ ( 1 , ⋯ , ), i.e., The integration operation in Eq. (2) means that (j ) can be achieved by integrating (j 1 , … , j ) over a n-dimensional hyperplane = 1 + ⋯ + . The GFRFs are nature extensions of the linear frequency response function to the nonlinear systems, because (j 1 , … , j ) is equal to the linear transfer function when = 1. Being different from linear systems, nonlinear systems have more complicated input-output relationship due to the effects of frequency multiplication and intermodulation, giving rise to the difficulties in identifying the contribution of input to output. To address this issue, the concept of NOFRF was proposed by Lang and Billings [23], making the input-output relationship of nonlinear systems much more transparent and physically meaningful. For a stable time-invariant nonlinear system, the NOFRFs were defined by i n n i Eq. (6) gives a compacted expression for (j ) in a manner similar to the output frequency response of a linear system, where (j ) is regarded as the input spectrum and (j ) is taken as the linear transfer function. Because (j ) can be formulated as the Fourier transform of the time sequence ( ), that is, ( ) can be seen as the input signals of the linear system (j ) in time domain. Substituting Eq. (6) into Eq. (2) gives which reflects that the nonlinear system under study can be decomposed into a set of linear systems. Therefore, the input-output relationship of the nonlinear system can be illustrated by a multiple-input/single-output (MISO) linear system as plotted in Figure 1. { ( )} are the inputs and { (j )} represent the linear transfer functions of the equivalent linear system. It should be noted that { (j )} depend not only on the nonlinearity of the system but also on the inputs as expressed in Eq. (4).

Figure 1 The MISO model of the nonlinear system linearized by NOFRFs
The frequency range of the nth order output in Figure 1 is determined by the hyperplane 12 n where ∈ [− , − ]⋃[ , ] , = 1, ⋯ , , 0 ≤ < , and [ , ] represents the non-negative frequency range of the original input ( ). In order to identify { (j )} and evaluate the contributions of each order nonlinearity to the output, a digital approach is proposed in the next section based on the CSA technique of MISO linear systems.

Application of CSA in Evaluating NOFRFs
As for the MISO linear system, the input-output relationship can be quantified by the ordinary coherence function under the condition that the inputs are uncorrelated, where (j ), (j ) and (j ) represent the auto-power and cross-power spectra and of the ith input and the output respectively. In the study, however, the inputs of the equivalent MISO linear system may be correlated with each other, resulting in the difficulties in distinguishing the effects of different order nonlinearities to the output. Conditioned spectral analysis provide a feasible way for solving the correlation problem of MISO linear systems, but it cannot be applied directly to nonlinear systems because the multiplication and intermodulation components of frequency are not taken into account in the existing CSA theories. Based on the NOFRFs model, the CSA technique will be extended to nonlinear systems in the study. As shown in Figure 1, a stable time-invariant nonlinear system with a general input can be decomposed into a MISO linear system near an equilibrium point. According to the procedure of CSA, the correlations between the inputs of different orders can be removed by taking conditioning operations, obtaining an equivalent MISO linear system with multiple uncorrelated inputs as shown in Figure 2. { •( −1)! (j ), = 1, ⋯ , + 1} represent the conditioned inputs, that is, the ith order input is conditioned by all of the lower order inputs { 1 (j ), 2 (j ), ⋯ , −1 (j )} in sequence, meaning that the parts of (j ) correlated with the lower order inputs are all removed. Thus, the conditioned inputs are uncorrelated with each other. The first order input 1 (j ) can be regarded as 1•0! (j ) where 0! denotes no conditions are operated. The last conditioned input ( +1)• ! (j ) refers to the measured output conditioned by the inputs { (j ), = 1, 2, ⋯ , } , equaling to the spectrum of the noise ( ) in Figure 1.
It is worth noting that the inputs are conditioned here from the first order to highest order, that is, the first input in Figure 2 is identical with the first one in Figure 1, and the th ( > 1) conditioned input equals to the th order input in Figure 1 conditioned by the lower orders. Obviously, conditioning in this sequence conforms to the physical meaning of Volterra series model of the nonlinear system, because the th input in Figure 1 depends on the order of frequency intermodulation and the higher order intermodulation is based on the intermodulation of lower orders.

Figure 2
The CSA model of the nonlinear system linearized by NOFRFs.
The conditioned inputs { •( −1)! (j ), 2 ≤ ≤ + 1} are formulated as which can be implemented by the iteration where 1 ≤ ≤ − 1 . The conditioned linear transfer functions { (j )} are given by in which 1 ≤ < ≤ + 1 , ( +1) (j ) refers to (j ), •( −1)! (j ) represents the auto-power spectral density of the conditioned input •( −1)! (j ) , and •( −1)! (j ) the cross-power spectral density of the conditioned inputs •( −1)! (j ) and •( −1)! (j ). The power spectral densities are defined by and respectively, where E[•] denotes the value of expectation and the superscript * signifies complex conjugation. The conditioned input power spectra can be obtained by and in the case of = , Iterating Eqs. (15), (18) and (19) (20) provides an alternative way for evaluating the NOFRFs of the nonlinear system. To sum up, the proposed digital method of evaluating NOFRFs of the nonlinear systems with a general input is carried out by the algorithm as shown in Figure 3. As a result, the output spectra generated by each order unconditioned input can be yielded by Eq. (6), such that the contributions of each order nonlinearity to the output are distinguished by this way. In addition, the total output power spectrum of the nonlinear system responding to the original input ( ) is given by (21) and the contribution ratio of the original input to the measured output can be quantified by

Simulation studies
In order to verify the effectiveness of the proposed method, two kinds of nonlinear systems are simulated in this section.

The first nonlinear system
The first nonlinear system is a known mechanical oscillator as shown in where ( ) is the force acted on the mass m， 0 ( ) is the displacement of the mass, ( ) represents the measured displacement of the mass, ( ) denotes a zero-mean Gaussian noise imposed on the measured data in which −2.56 ≤ ≤ 2.56 seconds as plotted in Figure  5(a). The output of the system is numerically simulated by using the fourth order Runge-Kutta method, as shown in Figure 5(b), and the input and output amplitude spectra obtained by FFT are sketched in Figure 6. The output spectrum distributes in the lower range of 0～20Hz besides the input frequency range of 30～55Hz, particularly near 15Hz where a distinct peak (i.e., the peak A) appears, implying that a part of power is transferred to the lower frequency range. The amplitude spectra of the system input and output obtained by FFT: (a) the input, (b) the output.

Correlations of the input between different orders
Given the time sequence of the input, the higher order input spectra { (j )} in Eq. (9) can be yielded by taking Fourier transform of the power series ( ). For example, the amplitude spectra of the first five order inputs obtained by FFT, i.e. {| (j ), = 1, ⋯ ,5|} , are as plotted in Figure 7. The magnitude and distribution of the five spectra are different obviously. Both the number and magnitude of crests increase with the order of input, and the input frequency range also widens with orders. is applied, where (j ) represents the cross-spectral density functions of the ith and jth order inputs, (j ) and (j ) represent the auto-power spectral density functions respectively. The results are as plotted in Figure  8, which indicate that there exist significant correlations between the inputs of different orders. Take { 1 2 , = 2, 3, 4, 5} for example, as shown in Figure 8(a), 1 ( ) and 2 ( ) are highly correlated in the two narrow bands, i.e. 24～30Hz and 55～62Hz, while 1 ( ) and 3 ( ) have correlation in the lower frequency range of 0～13Hz and the higher frequency range of 70～100Hz, and 1 ( ) and 4 ( ) have notable correlation in the middle range of 30 ～ 55Hz where the correlation between 1 ( ) and 5 ( ) is even higher. Besides, the correlation between 1 ( ) and 5 ( ) is remarkable in the ranges of 0～10Hz and 77～96Hz. Similarly, distinct correlations also exist between the 2nd, 3rd and 4th order inputs, which can be seen in Figure 8 The coherence spectra of the input between different orders: (a) between the first and the higher orders, (b) between the second and the higher orders, (c) between the third and the higher orders, (d) between the fourth and the fifth orders.

Nonlinear output response analysis based on CSA method
To reveal the contributions of each order nonlinearity to the system output, CSA technique is applied here. The conditioned input amplitude spectra {| •( −1)! (j )|, = 1, 2, ⋯ , 5} are shown in Figure 9, between which the correlations are all removed. The output power spectral components corresponding to the conditioned inputs are plotted in Figure 10 is displayed in Figure 10(f). It is illustrated that the output component correlated with the first order input is prominent in the frequency range of 30～55Hz, while the output component correlated with the second order conditioned input is dominated in the lower frequency range of 0～20Hz. Nevertheless, it cannot be stated by now that the output power distributed in 30～ 55Hz is generated mainly by the first order input, because all of the correlated output components induced actually by the second to fifth order nonlinearities of the system are attributed to the first order in the CSA procedure. Similarly, it cannot be concluded until now that the transferred power in 0～20Hz is generated mainly by the second order conditioned input, because all of the correlated output components induced actually by the third to fifth order nonlinearity are also attributed to the second order in the lower frequency range in the method.  The first five order NOFRFs { (j ), = 1, 2, ⋯ , 5} of the nonlinear system are obtained by Eq. (24), and their amplitude spectra {| (j )|} are plotted in Figure 11, which reflect the amplitude-frequency response properties of each order nonlinearity respecting the specific input. The second to fifth order NOFRFs all have higher amplitudes in the frequency range of 30～55 Hz, meaning that the nonlinear system has great gains on the amplitudes of the second to fifth order unconditioned inputs in the frequency range. Besides, both | 3 (j )| and | 5 (j )| exhibit many distinct peaks in the lower frequency range of 0～20 Hz, implying that the system also have notable gains on the amplitudes of the third and fifth order unconditioned inputs in the lower range. Nevertheless, because the ith order output frequency response is equal to the product of (j ) and (j ) as expressed by Eq. (6), the nonlinear output frequency responses of each order depends not merely on the NOFRF but also on the input. By Eq. (6) one can obtain exactly the output spectra of the nonlinear system responding to each order unconditioned input, which reveal the contributions of each order nonlinearity to the system output without considering interferences. The amplitude spectra of the unconditioned output responses, i.e. {| (j )|, = 1, ⋯ , 5}, are plotted in Figure 12(a)～(e). The first, third and fifth order inputs induce remarkable output responses in the frequency range of 30～55Hz which is identical with the first order input frequency range, while the second and fourth order inputs mainly give rise to the output components in the lower frequency range and the higher frequency range.
The amplitude spectrum of total output yielded by summing the unconditioned output response components, i.e. |∑ (j )

=1
|, is shown in Figure 12(f), which matches the result of Runge-Kutta method in Figure 6(b) very well, verifying numerically that the CSA based method is effective in distinguishing the contributions of each order nonlinearity to the output of the system. For example, at the frequency of 15.64Hz where the distinct peak A takes place in the range of 10 ～ 20Hz, the first five order responses and their summation are listed in Table 1. The magnitude of the first order output response at the frequency is only 5.4055 × 10 −9 m which is much smaller than other orders, so its influence on the output is slight. Although the magnitudes of the second and fourth order responses are higher at the frequency, but their phases are inverse, leading to the reduction amplitude. For the same reason, the output responses of the third and fifth orders are also cut down. As a result, the magnitude of the summation of the first five order output responses is equal to 1.1299 × 10 −5 m which is close to the numerical result 1.0993 × 10 −5 m obtained by Runge-Kutta method. The deviation of 2.784% is mainly due to the noise ( ) imposed on the output, which is common in practical systems.

The second nonlinear system
The second case for our simulation is a quadratic system with time lags as shown in Figure 13. Its governing equation is given by where ( ) and 0 ( ) represent the input and output of the system in time domain respectively, ( ) is the measured output of the system, ( ) denotes a zero-mean Gaussian noise acted on the output, 1 = −0.64, 2 = 1.0 , 3 = 0.9 and 4 = 1.0 are characteristic parameters, 1 = 0.001s and 2 = 0.0005s refer to the time delays. Suppose that both input and output of the system are displacements and the input is formed by Eq. (24) as same as the first nonlinear system. The output sequence of time obtained directly by Eq. (25) is shown in Figure 14(a), and the amplitude spectrum of the output obtained by FFT is sketched in Figure 14(b). Obviously, the magnitude of the output spectrum is prominent in the two ranges of 0～25Hz and 60～110Hz, while in the frequency range of input, i.e., 30～55Hz, the magnitude is approximately equal to zero, implying that almost all of the power is transferred to the lower and higher frequency ranges. Using the CSA method, the first five order NOFRFs { (j ), = 1, ⋯ , 5} of the second system are obtained and their amplitude spectra {| (j )|} are as shown in Figure 15, which reflect the amplitude-frequency response characteristics of each order nonlinearity of the system respecting to the input. Like the first nonlinear system simulated above, all the amplitude spectra of the second to fifth order NOFRFs have higher magnitudes in the frequency range of 30～55 Hz which is approximately equal to the original input frequency range. Besides, many distinct peaks exist in the lower and higher frequency ranges, e.g. | 3 (j )| in the range of 0～20 Hz, indicating that the NOFRFs also have notable gains on the input at the frequencies. However, because the ith order output frequency response is equal to the product of (j ) and (j ) as given by Eq. (6), the output frequency responses magnitude reduces to 41.4629 μm as a result of interferences. In addition, it is uncovered that the contribution of the first order is much smaller than the second, fourth and fifth orders at the frequency, and the contribution of the third order is even slighter. The total output magnitude at the frequency is approximately equal to the numerical result 40.7236μm obtained by the fourth order Runge-Kutta method with a deviation of 1.815% which is mainly due to the noise ( ) acting on the output.

Conclusions
(1) In this study, a digital method of evaluating NOFRFs of the nonlinear systems with a general input is proposed by applying CSA technique, i.e., the NOFRFs-CSA method as shown in Figure 3. (2) The interferences among various order output frequency components are evaluated and removed in the analysis procedure, and the NOFRFs are identified consequently. (3) Based on the NOFRFs-CSA method, the contributions of different order nonlinearities of the nonlinear system to the output can be distinguished in frequency domain, revealing the nonlinear frequency response characteristics of the nonlinear system. (4) Two kinds of nonlinear systems expressed by Eq. (23) and Eq. (25) were simulated with a same general input respectively. The results reached by the proposed method are very close to the numerical results obtained by the fourth order Runge-Kutta method, verifying the effectiveness of the proposed method on the evaluation of NOFRFs.