A Comparative Study of Parameter Identification Methods for Asymmetric Nonlinear Systems with Quadratic and Cubic Stiffness

Understanding the nonlinear dynamic characteristics of engineering structures is challenging, especially for the systems that exhibit asymmetric nonlinear behavior. This paper compared four parameter identification methods for asymmetric nonlinear systems incorporating quadratic and cubic stiffness nonlinearities. Hilbert transform, zero-crossing, direct quadrature, and wavelet transform were used to obtain the backbone, envelope, and restoring force curves from the free vibration time history. A nonlinear curve-fitting method was then applied to estimate the stiffness parameters of the asymmetric systems, and a linear least square fitting approach was utilized to estimate the damping parameters of the asymmetric systems. We used the Helmholtz–Duffing oscillator as a numerical example and a nonlinear vibration absorber with geometric imperfections to verify the feasibility and accuracy of these methods. The advantages and disadvantages of these methods and the deviations in estimated results are discussed.


Introduction
Asymmetric nonlinear engineering structures, such as a mistuning quasi-zero-stiffness vibration isolator [1,2], cables [3], and geometric imperfect plate-like designs [4][5][6][7], have attracted widespread attention. The dynamic characteristics of these asymmetric systems are more complicated to analyze than those of the symmetric systems. For example, constant drift can occur in the response when there is little linear stiffness in the quasizero-stiffness vibration isolator. The nonlinear isolation system exhibits a mixed softening and hardening characteristic [1], which results in the multiple jump phenomena and hysteretic behavior [2]. Multivaluedness of the response curves occurs with different features depending on cables (or plates) and excitation force parameters [3][4][5][6][7]. To fully understand the nonlinear aspect of asymmetric systems, nonlinear parameter identification is one of the crucial procedures [8][9][10].
In recent years, backbone curves have been used to identify stiffness parameters of asymmetric systems. The calculation methods of the backbone curves can be classified into two types: analytical and numerical methods. Analytical methods include harmonic balance, multi-scale [11], and normal form methods [12,13]. The comparison of these analytical methods can be referred to [14]. Some software packages have implemented numerical algorithms based on a nonlinear normal mode framework [15][16][17]. Common experimental methods for extracting backbone curves include the resonance decay method [18], the control-based continuation method [19,20] and phase-locked loops [21]. For multi-degreeof-freedom systems, the nonlinear normal modes of interest are usually isolated by the force appropriation method [22][23][24]. The resonance decay method is then used to estimate m .. (1) where m is the mass, c is the damping coefficient, f (t) is the excitation force, k 1 , k 2 , and k 3 are the linear, quadratic, and cubic nonlinear stiffness coefficients, respectively. If the system vibrates freely with an initial displacement, f (t) = 0. Then, dividing both sides of Equation (1) by the mass m yields ..
x + ω 2 0 (t)x = 0 (2) where h = c/2m is the damping factor and ω 0 (t) = (k 1 + k 2 x + k 3 x 2 )/m is the instantaneous frequency. They have the units of 1/s. It can be seen that the analytical solution for the viscous damping force per unit mass is written as f c (x) = 2h .
x when the range of velocity is given, and the analytical solution for the restoring force per unit mass is given by where the regime of displacement is determined. Backbone curves for the Helmholtz-Duffing oscillator are more complex than those for the symmetric systems [42]. Using the harmonic balance method and substituting the approximate harmonic solution x(t) = A 0 + A 1 cos(ωt + θ) into Equation (1) yield where 3 is the bias term, A 2 1 = − k 3 A 3 0 + k 2 A 2 0 + k 1 A 0 /(3k 3 A 0 /2 + k 2 /2) is the first harmonic term, β = k 1 − k 2 2 /3k 3 , and f 0 = k 1 k 2 /3k 3 − 2k 3 2 /27k 2 3 . Adding more harmonic terms to the approximate harmonic solution can yield a more accurate solution, but the calculation of the backbone curve is complicated. Substituting the approximate harmonic solution x(t) = A 0 + A 1 cos(ωt + θ) + A 2 cos(2ωt + ϕ), where the second harmonic term A 2 is included, into Equation (1) for the lightly damped case, gives which can be compared well with the numerical results of the backbone curves [42]. The comparison between analytical solutions and the numerical results is shown in Section 3.

Identification Methods for the Characteristic Curves
This section introduces four methods for obtaining the characteristic curves of the Helmholtz-Duffing oscillator, which are restoring force curves, damping force curves, envelopes, and backbone curves. These curves are combined with the analytical solutions given in Section 2 to estimate the stiffness and damping parameters of asymmetric structures. The details are discussed in Sections 4 and 5. For the sake of simplicity, the parameters for the oscillator are m = 0.1 kg, c = 0.4 Ns/m, k 1 = 4000 N/m, k 2 = −10 7 N/m 2 , and k 3 = 10 10 N/m 3 . The initial displacement and velocity are x(0) = 0.0018 m and .
x(0) = 0 m/s. The free-decay response is numerical integration calculated by using the fourth-order Runge-Kutta method, and the sampling frequency is f s = 2000 Hz. For the asymmetric system, when the free decay response is measured, the restoring force curve is constructed separately from the positive and negative signal parts [35], which is given by where A c and ω c included positive and negative congruent envelopes and congruent natural frequency. The viscous damping force is approximately given by where h(t) is the instantaneous damping factor and A .
x (t) is the congruent envelope of velocity. The backbone curves can be obtained using the relationship between instantaneous frequency and amplitude of each harmonic term.

Hilbert Transform
Hilbert transform (HT) has been widely used in the parameter identification and signal decomposition of nonlinear systems. Feldman applied this method to identify free and forced vibration systems, and further applied the nonparametric identification method to the asymmetric systems [43][44][45][46]. The complex analytic form of a free-decay response is given by X(t) = x(t) + j x(t) = A(t)e jϕ(t) , where x(t) is the Hilbert transform of the signal x(t), the envelope A(t) = x 2 (t) + x 2 (t), and the instantaneous phase . The instantaneous undamped natural frequency is given by and the instantaneous damping factor is given by ω are the first and second derivatives of envelope and frequency, respectively. Hilbert vibration decomposition (HVD) is a time-varying vibration decomposition method based on the Hilbert transform [45]. The main harmonic components of asymmetric systems can be obtained by using the HVD. By weighted summing the decomposed harmonic components, the congruent envelope, and the congruent frequency can be obtained. The congruent envelope A c is given by where A l is the envelope of the lth order component, and φ l is the phase angle between the primary component and lth order component. Instantaneous natural frequency ω 0 can be decomposed into a sum of high-order intrinsic components. The congruent natural frequency ω c is given by where ω 0l is the envelope of the lth order instantaneous natural frequency, and φ ωl is the phase angle between the primary component and lth order component. The Hilbert transform method and the relevant Matlab programs [46] are used in this paper. Free-decay response of the Helmholtz-Duffing oscillator and its envelopes are shown in Figure 1a, where the positive and negative congruent envelopes are obtained by Equation (10). The instantaneous natural frequency and their envelopes are obtained by Equations (8) and (11), as shown in Figure 1b. Figure 1c shows the first fourth components of the free vibration obtained by using HVD. In order to remove the end effect of the Hilbert transform, only the analyzed results between 0.2 to 1.7 s are chosen. The restoring force curves are shown in Figure 1d, where the analytical solution is given by Equation (3), and the numerical results are obtained by Equation (6). The backbone curves of the first harmonic term A 1 and bias term A 0 are plotted in Figure 1e-f, respectively, where the analytical solutions are given by Equation (5), and the numerical results are obtained by using HVD. The logarithmic form of the envelope, instantaneous damping factor obtained by Equation (9), and damping force curve obtained by Equation (7) are shown in Figure 1g-i, respectively. For a weakly nonlinear system, the analytical solution of the logarithmic envelope is approximately given by −ct/2m + ln x(0), as shown in Figure 1g. The damping factor is h = c/2m = 21/s and the analytical damping force is f c (x) = 2h .
x, which are plotted in Figure 1h,i, respectively.

Zero-Crossing
The zero-crossing (ZC) method is a simple method for estimating the instantaneous frequency and amplitude of free decay signals [27,36]. Londoño et al. [27] combined the zero-crossing and peak picking methods to obtain the backbone curves of symmetric nonlinear systems. Ondra et al. [36] extended the zero-crossing method to obtain the backbone and restoring force curves of asymmetric systems. This method is explained in detail in Figure 2. The ith zero-crossing point t i , the pth positive peak points A p , and the nth negative peak points A n are first obtained using the zero and peak picking procedure. The instantaneous frequency at the zero-crossing point t i is The positive and negative envelopes at the zero-crossing points are then obtained by linear interpolating the positive and negative peak points, respectively. The corresponding instantaneous frequencies for positive and negative envelope points are ω p = π/T p and ω n = π/T n , respectively. using HVD. The logarithmic form of the envelope, instantaneous damping factor obtained by Equation (9), and damping force curve obtained by Equation (7) are shown in Figure  1g-i, respectively. For a weakly nonlinear system, the analytical solution of the logarithmic envelope is approximately given by −+ 2 ln (0) ct m x , as shown in Figure 1g. The damping factor is h = c/2m = 2 1/s and the analytical damping force is = ( ) 2 c f x hx , which are plotted in Figure 1h,i, respectively.

Zero-Crossing
The zero-crossing (ZC) method is a simple method for estimating the instantaneous frequency and amplitude of free decay signals [27,36]. Londoño et al. [27] combined the zero-crossing and peak picking methods to obtain the backbone curves of symmetric nonlinear systems. Ondra et al. [36] extended the zero-crossing method to obtain the backbone ing instantaneous frequencies for positive and negative envelope points are ωp = π/Tp and ωn = π/Tn, respectively.
The amplitudes of the first harmonic term A1 and bias term A0 are given by The free decay response is given in Figure 3a again. Then, the zero-crossing points and positive and negative peak points can be obtained and plotted in Figure 3a. The instantaneous frequencies for the positive and negative envelopes are shown in Figure 3b. Substituting the positive and negative envelopes and instantaneous frequencies in Equation (6) yields the restoring force curves, as shown in Figure 3c. It should be noticed that the time history response of asymmetric systems normally contains multiple harmonic components. Unlike the HVD, the backbone curves in the zero-crossing method are calculated using Equations (12) and (13), where the high order harmonic terms induce estimation errors, especially for the bias term. Before calculating the instantaneous frequency and amplitudes, the high-order harmonic terms should be filtered to make sure that the response can be approximately written as For the numerical example in this section, the free decay signal is passed through a low-pass filter with a cutoff frequency of 50 Hz. The backbone curves are shown in Figure 3d,e. It can be seen that the backbone curves obtained from the filtered response are closer to the analytical backbone curves. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 3f-h. Numerical results for the damping factor and damping force are obtained by Equations (7) and (9), where the envelope of velocity is obtained by the peak picking method. To unify several methods, the data used in the ing instantaneous frequencies for positive and negative envelope points are ωp = π/Tp and ωn = π/Tn, respectively. The amplitudes of the first harmonic term A1 and bias term A0 are given by The free decay response is given in Figure 3a again. Then, the zero-crossing points and positive and negative peak points can be obtained and plotted in Figure 3a. The instantaneous frequencies for the positive and negative envelopes are shown in Figure 3b. Substituting the positive and negative envelopes and instantaneous frequencies in Equation (6) yields the restoring force curves, as shown in Figure 3c. It should be noticed that the time history response of asymmetric systems normally contains multiple harmonic components. Unlike the HVD, the backbone curves in the zero-crossing method are calculated using Equations (12) and (13), where the high order harmonic terms induce estimation errors, especially for the bias term. Before calculating the instantaneous frequency and amplitudes, the high-order harmonic terms should be filtered to make sure that the response can be approximately written as For the numerical example in this section, the free decay signal is passed through a low-pass filter with a cutoff frequency of 50 Hz. The backbone curves are shown in Figure 3d,e. It can be seen that the backbone curves obtained from the filtered response are closer to the analytical backbone curves. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 3f-h. Numerical results for the damping factor and damping force are obtained by Equations (7) and (9), where the envelope of velocity is obtained by the peak picking method. To unify several methods, the data used in the , zero-crossing points; ing instantaneous frequencies for positive and negative envelope points are ωp = π/Tp and ωn = π/Tn, respectively. The amplitudes of the first harmonic term A1 and bias term A0 are given by The free decay response is given in Figure 3a again. Then, the zero-crossing points and positive and negative peak points can be obtained and plotted in Figure 3a. The instantaneous frequencies for the positive and negative envelopes are shown in Figure 3b. Substituting the positive and negative envelopes and instantaneous frequencies in Equation (6) yields the restoring force curves, as shown in Figure 3c. It should be noticed that the time history response of asymmetric systems normally contains multiple harmonic components. Unlike the HVD, the backbone curves in the zero-crossing method are calculated using Equations (12) and (13), where the high order harmonic terms induce estimation errors, especially for the bias term. Before calculating the instantaneous frequency and amplitudes, the high-order harmonic terms should be filtered to make sure that the response can be approximately written as For the numerical example in this section, the free decay signal is passed through a low-pass filter with a cutoff frequency of 50 Hz. The backbone curves are shown in Figure 3d,e. It can be seen that the backbone curves obtained from the filtered response are closer to the analytical backbone curves. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 3f-h. Numerical results for the damping factor and damping force are obtained by Equations (7) and (9), where the envelope of velocity is obtained by the peak picking method. To unify several methods, the data used in the , positive peak points; ing instantaneous frequencies for positive and negative envelope points are ωp = π/Tp and ωn = π/Tn, respectively. The amplitudes of the first harmonic term A1 and bias term A0 are given by The free decay response is given in Figure 3a again. Then, the zero-crossing points and positive and negative peak points can be obtained and plotted in Figure 3a. The instantaneous frequencies for the positive and negative envelopes are shown in Figure 3b. Substituting the positive and negative envelopes and instantaneous frequencies in Equation (6) yields the restoring force curves, as shown in Figure 3c. It should be noticed that the time history response of asymmetric systems normally contains multiple harmonic components. Unlike the HVD, the backbone curves in the zero-crossing method are calculated using Equations (12) and (13), where the high order harmonic terms induce estimation errors, especially for the bias term. Before calculating the instantaneous frequency and amplitudes, the high-order harmonic terms should be filtered to make sure that the response can be approximately written as For the numerical example in this section, the free decay signal is passed through a low-pass filter with a cutoff frequency of 50 Hz. The backbone curves are shown in Figure 3d,e. It can be seen that the backbone curves obtained from the filtered response are closer to the analytical backbone curves. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 3f-h. Numerical results for the damping factor and damping force are obtained by Equations (7) and (9), where the envelope of velocity is obtained by the peak picking method. To unify several methods, the data used in the , negative peak points; ing instantaneous frequencies for positive and negative envelope points are ωp = π/Tp and ωn = π/Tn, respectively. The amplitudes of the first harmonic term A1 and bias term A0 are given by The free decay response is given in Figure 3a again. Then, the zero-crossing points and positive and negative peak points can be obtained and plotted in Figure 3a. The instantaneous frequencies for the positive and negative envelopes are shown in Figure 3b. Substituting the positive and negative envelopes and instantaneous frequencies in Equation (6) yields the restoring force curves, as shown in Figure 3c. It should be noticed that the time history response of asymmetric systems normally contains multiple harmonic components. Unlike the HVD, the backbone curves in the zero-crossing method are calculated using Equations (12) and (13), where the high order harmonic terms induce estimation errors, especially for the bias term. Before calculating the instantaneous frequency and amplitudes, the high-order harmonic terms should be filtered to make sure that the response can be approximately written as For the numerical example in this section, the free decay signal is passed through a low-pass filter with a cutoff frequency of 50 Hz. The backbone curves are shown in Figure 3d,e. It can be seen that the backbone curves obtained from the filtered response are closer to the analytical backbone curves. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 3f-h. Numerical results for the damping factor and damping force are obtained by Equations (7) and (9), where the envelope of velocity is obtained by the peak picking method. To unify several methods, the data used in the , envelopes at zero-crossing points obtained by interpolation.
The amplitudes of the first harmonic term A 1 and bias term A 0 are given by The free decay response is given in Figure 3a again. Then, the zero-crossing points and positive and negative peak points can be obtained and plotted in Figure 3a. The instantaneous frequencies for the positive and negative envelopes are shown in Figure 3b. Substituting the positive and negative envelopes and instantaneous frequencies in Equation (6) yields the restoring force curves, as shown in Figure 3c. It should be noticed that the time history response of asymmetric systems normally contains multiple harmonic components. Unlike the HVD, the backbone curves in the zero-crossing method are calculated using Equations (12) and (13), where the high order harmonic terms induce estimation errors, especially for the bias term. Before calculating the instantaneous frequency and amplitudes, the high-order harmonic terms should be filtered to make sure that the response can be approximately written as x(t) ≈ A 0 + A 1 cos(ωt + θ). For the numerical example in this section, the free decay signal is passed through a low-pass filter with a cut-off frequency of 50 Hz. The backbone curves are shown in Figure 3d,e. It can be seen that the backbone curves obtained from the filtered response are closer to the analytical backbone curves. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 3f-h. Numerical results for the damping factor and damping force are obtained by Equations (7) and (9), where the envelope of velocity is obtained by the peak picking method. To unify several methods, the data used in the restoring force curve, logarithmic envelope, and damping curve are also approximately taken from 0.2 to 1.7 s. The data used in backbone curves are relatively short because part of the response in the large amplitude regime is filtered out. The analytical solutions shown in Figure 3 are obtained using the similar approaches mentioned in Section 3.1.

Direct Quadrature
The direct quadrature (DQ) method was proposed by Huang et al. [38]. Firstly, because the Hilbert transform of a product of functions is limited by the Bedrosian theorem [47], a normalization scheme was proposed to separate amplitude modulation (AM) and frequency modulation (FM) of the signal. Secondly, according to the Nuttall theorem [48], it is not applicable for all signals to obtain their quadrature forms using the Hilbert transform. Therefore, the direct quadrature method is used. The direct quadrature method has been applied to the symmetric signal. However, when it is applied to the asymmetric signal, some modifications are required. restoring force curve, logarithmic envelope, and damping curve are also approximately taken from 0.2 to 1.7s. The data used in backbone curves are relatively short because part of the response in the large amplitude regime is filtered out. The analytical solutions shown in Figure 3 are obtained using the similar approaches mentioned in sub-Section 3.1.

Direct Quadrature
The direct quadrature (DQ) method was proposed by Huang et al. [38]. Firstly, because the Hilbert transform of a product of functions is limited by the Bedrosian theorem [47], a normalization scheme was proposed to separate amplitude modulation (AM) and frequency modulation (FM) of the signal. Secondly, according to the Nuttall theorem [48], it is not applicable for all signals to obtain their quadrature forms using the Hilbert transform. Therefore, the direct quadrature method is used. The direct quadrature method has The first step is to find the positive and negative peak points and use the cubic spline function to obtain the positive and negative envelopes A p and A n . Then the normalization process is carried out. The positive time-domain response is normalized by y p (t p ) = x p (t p )/A p and the negative time-domain response is normalized by y n (t n ) = x n (t n )/A n , where x p (t p ) and x n (t n ) are the positive and negative time-domain responses, y p (t p ) and y n (t n ) are the normalized positive and negative responses. Repeat the above steps until the normalized results are all in [−1,1]. The FM part of the signal is F(t) = y l (t) = y l p (t p ), y l n (t n ) , where l is the number of iterations. After normalization, the AM part is given by , A n (t n ) = x n (t n ) y l n (t n ) (14) Then using the cubic spline function, the positive envelope A p (t) and negative envelope A n (t) in the entire time domain are obtained. The amplitudes of the first harmonic and the bias terms are given by The time domain response and its positive and negative envelopes are shown in Figure 4a. The FM part can be regarded as a sinusoid, so the instantaneous frequency is The FM part whose absolute values are less than 0.9 is used to calculate the instantaneous frequency using Equation (16), and the rest of the instantaneous frequency points are interpolated using a cubic spline. The result is shown in Figure 4b. The restoring force curve obtained by Equation (6) is shown in Figure 4c. Similar to the zero-crossing method, the amplitudes of the first harmonic and the bias terms are also sensitive to the high-frequency components of the free decay response. Before calculating the backbone curves, the free decay response is passed through a low-pass filter with a cut-off frequency of 50 Hz. The estimated backbone curves of A 1 and A 0 are shown in Figure 4d,e, where amplitudes are obtained by Equation (15), and the frequency is the filtered instantaneous frequency. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 4f-h. The analytical solutions are also shown in Figure 4.

Wavelet Transform
Wavelet transform (WT) is a time-frequency analysis tool that can automatically adjust the size of the analysis window with the change of frequency. With the development in recent years, wavelet analysis has been widely applied in nonlinear system identification [39,[49][50][51].
For the free decay response of the Helmholtz-Duffing oscillator shown in Figure 5a, the frequency spectrum can be obtained from the Matlab cwt function using the Morlet wavelet and is shown in Figure 5b. From this figure, the envelope of the signal is given by where a is the scale parameter, b is the translation parameter, |W x (a(b), b)| = max a |W x (a, b)| is the maximum value of wavelet coefficients at each time point. The instantaneous frequency is obtained from the frequency points corresponding to the maximum value of the wavelet coefficients at each time point. In order to obtain the smooth envelope and instantaneous frequency, the results are passed through a low-pass filter with a cut-off frequency of 20 Hz. The results are shown in Figure 5a,b. The backbone curve of A 1 constructed by the instantaneous amplitude and frequency is shown in Figure 5c. The bias term is obtained by the wavelet decomposition using the Matlab function wavedec, as shown in Figure 5d. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 5e-g. It can be seen that the numerical solutions can be reasonably compared well with the analytical solutions.

Wavelet Transform
Wavelet transform (WT) is a time-frequency analysis tool that can automatically adjust the size of the analysis window with the change of frequency. With the development in recent years, wavelet analysis has been widely applied in nonlinear system identification [39,[49][50][51].
For the free decay response of the Helmholtz-Duffing oscillator shown in Figure 5a, the frequency spectrum can be obtained from the Matlab cwt function using the Morlet wavelet and is shown in Figure 5b. From this figure, the envelope of the signal is given by where a is the scale parameter, b is the translation parameter, is the maximum value of wavelet coefficients at each time point. The instantaneous frequency is obtained from the frequency points corresponding to the maximum value of the wavelet coefficients at each time point. In order to obtain the smooth envelope and instantaneous frequency, the results are passed through a low-pass filter with a cut-off frequency of 20 Hz. The results are shown in Figure 5a,b. The backbone curve of A1 constructed by the instantaneous amplitude and frequency is shown in Figure 5c. The bias term is obtained by the wavelet decomposition using the Matlab function wavedec, as shown in Figure 5d. The logarithmic form of the envelope, instantaneous damping factor, and damping force curve are shown in Figure 5e-g. It can be seen that the numerical solutions can be reasonably compared well with the analytical solutions.

Parameter Estimation and Discussion
In this section, the stiffness and damping parameters of the asymmetric systems are estimated from the characteristic curves of the Helmholtz-Duffing oscillator obtained in Section 3. The stiffness parameters are obtained by polynomial fitting the restoring force curve shown in Equation (3). The Matlab function polyfit computes the least square polynomial. This method is called the restoring force curve method (RFCM). Although the backbone curve of the first harmonic term obtained by using Equation (4) deviates from the numerical result obtained by using HVD in the bending regime [42], which is the curved part, we can also estimate the stiffness parameters from this curve. Matlab function fminsearch is used here to find the optimal stiffness parameters. This identification method is called the backbone curve method (BCM). The estimated stiffness parameters are shown in Table 1.
The viscous damping can be estimated by linear fitting the damping force curve, called the damping force curve method (DFCM). For a weakly nonlinear system, the envelope is approximately given by ln A(t) = −ct/2m + ln A 0 . Therefore, the natural logarithm of the envelope can also be used to estimate the damping [28], called the logarithmic envelope method (LEM). Matlab function polyfit is utilized for computing the least square linear coefficient. The estimated damping coefficients are shown in Table 2.  The estimated stiffness parameters are obtained by using the four methods. The Hilbert transform and Hilbert vibration decomposition combined with the low pass filter can give accurate estimation results. This method can not only decompose the signal and obtain the backbone curve of each harmonic term, but also combine harmonic terms to construct the restoring force. However, the Hilbert transform method has an end effect. The data at the beginning and end need to be removed. For zero-crossing with the peak picking method, it is straightforward to implement. Even for the asymmetric signal, this method can analyze the positive and negative time domain separately and construct the corresponding restoring force. However, when this method is applied to the signal contaminated with noise, the zero-crossing points and peak points are difficult to obtain accurately. The signal should be properly filtered to solve this issue. Also, there are only fewer points to extract for the short-time signal. In this case, the interpolation method can be used in the whole time domain to obtain more points. For the direct quadrature method, part of instantaneous frequency points are interpolated using a cubic spline, so the obtained instantaneous frequency is not accurate enough, especially for the signal with a low sampling rate or a large amount of normalized data is over 0.9. For the wavelet transform, the wavelet function needs to be selected carefully and appropriately.
The nonlinear stiffness parameters estimated by the restoring force curve using the zero-crossing and direct quadrature methods are not accurate. The estimated results of the backbone curve method are all well. The backbone curve method estimates the entire backbone curve, so the estimated stiffness parameters are comprehensively affected by the deviation of each amplitude regime.
For the damping coefficient, the four methods seem to achieve similar results. The estimated solutions of the logarithmic envelope method are all less than the actual value because it uses the analytical solution of the linear system. But this method is simple and easy to estimate for weakly damped systems. The error of the damping force curve method is small. The derivatives of the envelope, instantaneous frequency, and displacement should be obtained for the damping force curve method, so the result is disturbed by noise easily. A proper filter can be used to deal with the influence of the noise, and the Bayesian approach is a good way to measure the uncertainty of the identification results [29].

Experimental Description
The test rig is shown in Figure 6, where Figure 6a is the photo of the test rig, and Figure 6b,c are the elevation and plan views of the nonlinear vibration absorber. The vibration absorber consists of a 4.86 g mass attached to a thin circle brass plate of 0.15 mm thickness. Figure 6c also shows the contour plot for the measured geometric imperfections obtained by moving the laser sensor through the translation surface of the plate. It can be seen that the plate is not flat, and has a certain initial deflection, as shown in Figure 6b. method is small. The derivatives of the envelope, instantaneous frequency, and displacement should be obtained for the damping force curve method, so the result is disturbed by noise easily. A proper filter can be used to deal with the influence of the noise, and the Bayesian approach is a good way to measure the uncertainty of the identification results [29].

Experimental Description
The test rig is shown in Figure 6, where Figure 6a is the photo of the test rig, and Figure 6b,c are the elevation and plan views of the nonlinear vibration absorber. The vibration absorber consists of a 4.86 g mass attached to a thin circle brass plate of 0.15 mm thickness. Figure 6c also shows the contour plot for the measured geometric imperfections obtained by moving the laser sensor through the translation surface of the plate. It can be seen that the plate is not flat, and has a certain initial deflection, as shown in Figure 6b.  The excitation signal was generated by the LDS V406 shaker, then measured by the B&K 4517 accelerometer and Microtrak™ 3 LTS-050-10 laser sensor, respectively. The sampling frequency was 2000 Hz. Because the linear natural frequency of the vibration absorber was much higher than the natural frequency of the shaker, the absorber and the shaker can be regarded as a single-degree-of-freedom system [28]. The excitation signal supplied to the shaker was from Agilent 33512B signal generator and passed through an LDS PA500L power amplifier. The measured signals were sampled by a NI PXIe-4492 acquisition system after passing through a B&K 2693 conditioner. The equation of motion of the plate with geometric imperfections can be simplified as a Helmholtz-Duffing oscillator, as discussed in [4,28]. Therefore, the mathematical model of the free-decay response of the nonlinear vibration absorber is given by where m EQ = (m s + m v )(m a + m)/(m s + m v + m a + m) = 5.44 g is the equivalent mass of the system. m s , m v , m a , and m are the mass of armature, support structure of the absorber, accelerometer, and absorber mass, respectively. x is the vibration response of the experimental system. Damping c, linear stiffness k 1 , quadratic and cubic nonlinear stiffness k 2 and k 3 are the parameters to be estimated.

Estimation and Discussion
Before the free decay experiment, several slow frequency sweep experiments from low to high frequency were carried out to determine the system's jump-down frequency, which was about 201 Hz. Therefore, the excitation signal was switched off at 200 Hz. The circular plate exhibited large stiffness nonlinearity, and the vibration modes of the system were well-separated. In order to exclude the influence of the higher harmonics and high order modes, the displacement of the mass measured by the laser sensor was passed through a low pass filter with a cut-off frequency of 270 Hz, as shown in Figure 7a. The specific experimental instruments and testing procedures were described in [42]. The data used in identification is approximately taken from 0.205 to 0.37 s. Figure 7b shows the restoring force curve. Because the response decayed fast, cubic spline interpolation was used to construct envelope and frequency at each time point for the zero-crossing method. The first several large-amplitude points were removed due to the deviation. The restoring force was estimated by using Equation (3). The estimated stiffness parameters are shown in Table 3. The estimated results of the zero-crossing and direct quadrature methods are close to each other.  Figure 7c shows the backbone curve. The backbone curve obtained using the Hilbert transform is compared well with the backbone curve obtained using the zero-crossing method. The backbone curve obtained by using the wavelet transform shows a slight deviation from the above two methods, while the backbone curve obtained by using the direct quadrature method shows the most significant deviation. The estimated stiffness parameters are shown in Table 3. It can be seen that, except for the direct quadrature method, the estimated results of other methods are compared well. The identification results of the Hilbert transform can be used to reasonably reconstruct the experimental system's restoring force and backbone curves. The stiffness parameters are also estimated well by the zero crossing method and wavelet transform using the backbone curve.
The natural logarithm of the envelope and the damping force curve are shown in Figure 7d,e. The results of the four methods are similar, as shown in Table 4. It can be seen that the damping coefficient estimated by the logarithmic envelope method is less than that estimated by the damping force curve method. The logarithmic envelope method uses the analytical solution of the envelope of the linear system. However, the damping of the experimental system is not strictly linear. It can be noticed from the logarithmic envelope and damping force curve that the damping factor is amplitude dependant and decreases with time.  Figure 7b shows the restoring force curve. Because the response decayed fast, cubic spline interpolation was used to construct envelope and frequency at each time point for the zero-crossing method. The first several large-amplitude points were removed due to the deviation. The restoring force was estimated by using Equation (3). The estimated stiffness parameters are shown in Table 3. The estimated results of the zero-crossing and direct quadrature methods are close to each other. To determine whether the identification procedure is successful, a reconstructed free decay response is compared with the measured response as shown in Figure 7a. The approximate solution of Equation (18) is given by where A 0 (t) is the time-dependent amplitude (or envelope) of bias term, A 1 (t) is the time-dependent amplitude of the first harmonic term, and φ(t) is the time-dependent phase. The time-dependent phase is obtained by integrating the time-dependent damped natural frequency The time-dependent damped natural frequency, namely the backbone curve, is given by Equation (4). The parameters obtained by the Hilbert transform method using the backbone curve and the measured envelopes are used in reconstructing the free decay response. The results are shown in Figure 8a. For the small amplitude regime, using the multiple-scales method [11], the time domain response is approximately given by where 0 () At is the measured envelope of bias term, A1 is the initial amplitude of the first harmonic term.  = 1 EQ n km and  = EQ 1 2 c m k are the undamped natural frequency and the damping ratio of the underlying linear system respectively. Damping obtained by the logarithmic envelope method is used in reconstructing the time response.
The results are shown in Figure 8b. It can be seen that the reconstructed responses obtained by the two reconstruction methods are compared well with the measured responses. The reconstruction method is described in detail in [28], which is applied to identify a Duffing type nonlinear vibration absorber.

Conclusions
This paper compared four identification methods for the stiffness and damping parameters of asymmetric systems with square and cubic nonlinearities. We verified the feasibility and accuracy of these methods by a Helmholtz-Duffing numerical example and a nonlinear vibration absorber with geometric imperfections. Hilbert vibration decomposition decomposes the asymmetric signal to obtain the backbone curve of each harmonic term. The asymmetric restoring force curve is constructed by the positive and negative congruent envelopes and frequencies obtained by using a weighted summing of the decomposed harmonic components. The obtained curves compare well with the analytical solution, but the disadvantage of the Hilbert transform is the end effect. Zero-crossing For the small amplitude regime, using the multiple-scales method [11], the time domain response is approximately given by x(t) = A 0 (t) + A 1 e −ζω n t cos ω n t + 3 16ζ where A 0 (t) is the measured envelope of bias term, A 1 is the initial amplitude of the first harmonic term. ω n = k 1 /m EQ and ζ = c/2 m EQ k 1 are the undamped natural frequency and the damping ratio of the underlying linear system respectively. Damping obtained by the logarithmic envelope method is used in reconstructing the time response. The results are shown in Figure 8b. It can be seen that the reconstructed responses obtained by the two reconstruction methods are compared well with the measured responses. The reconstruction method is described in detail in [28], which is applied to identify a Duffing type nonlinear vibration absorber.

Conclusions
This paper compared four identification methods for the stiffness and damping parameters of asymmetric systems with square and cubic nonlinearities. We verified the feasibility and accuracy of these methods by a Helmholtz-Duffing numerical example and a nonlinear vibration absorber with geometric imperfections. Hilbert vibration decomposition decomposes the asymmetric signal to obtain the backbone curve of each harmonic term. The asymmetric restoring force curve is constructed by the positive and negative congruent envelopes and frequencies obtained by using a weighted summing of the decomposed harmonic components. The obtained curves compare well with the analytical solution, but the disadvantage of the Hilbert transform is the end effect. Zero-crossing with the peak picking method extracts positive and negative peak points and obtains instantaneous frequency more simply. The bias term and first harmonic term are obtained by using the sum and difference of the positive and negative envelopes. The direct quadrature method also analyzes the positive and negative signal parts separately and uses the same method as the zero-crossing method to obtain the bias term and the first harmonic term. The two methods are more sensitive to noise. It is essential for the wavelet transform to select wavelet function appropriately.
The nonlinear stiffness parameters estimated are not accurate in the numerical example when the restoring force curve is obtained using the zero-crossing method and direct quadrature method. The identification errors for nonlinear stiffness parameters are about 20%. The other methods can estimate the stiffness parameters accurately, and the identification errors are less than 5%. For the experimental system, the identification results of the Hilbert transform can describe the experimental system. The stiffness parameters estimated by the restoring force curve using the zero-crossing method and the direct quadrature method are close to each other. The backbone curve obtained by using the direct quadrature method deviates from those obtained by the other methods.
The logarithmic envelope can be used to estimate damping. The estimated results are the same, and all are less than the actual values. The other tool is the damping force curve, which is constructed by the damping factor and envelope of the velocity. The estimated results are better than those of the logarithmic envelope. However, the damping force curve needs the derivative of the envelope, instantaneous frequency, and displacement, so it has poor robustness against noise. It can be seen from the logarithmic envelope and damping force curve that the damping of the experimental system is not strictly linear.
This paper is to identify the parameters of a predetermined system model. Restoring the force curve method and damping force curve method are nonparametric identification methods. They produce the best functional representation of the system without a priori assumption about the system model. However, it is necessary for the backbone curve method to know the nonlinear system model and the theoretical solution of the backbone curve in advance. The logarithmic envelope method is suitable for lightly damped systems and is an approximate estimation method.