On the identification of damping from non-stationary free decay signals using modern signal processing techniques

A numerical implementation of system identification from non-linear and non-stationary signals is presented. The continuous wavelet transform (CWT) along with the complex Morlet wavelet skeleton curve extraction and Hilbert Transform (HT)-based methodologies are used for identification purposes. A comparison of the advantages of each technique in the analysis of non-stationary free decay systems is presented and improvements to the current methodologies are proposed. The HT approach offered good results in the estimation of the instantaneous amplitude in low damping and non-noisy signals. However, it is highly sensitive to impulses and irregularities in the signal, which affects the proper detection of frequency and amplitude parameters in real-life signals. The CWT exhibited better results for the analysis of noisy signals, from the resulting wavelet map the noise content can be distinguished from the actual system response. That is, the modes show a distinctive pattern in the map allowing proper modal extraction. However, for highly damped non-stationary decaying signals, the results are affected by the decay rate, round-up errors, and edge effects.


Introduction
In vibration-based structural health evaluation, damage is usually associated with changes or irregularities in the system dynamic response parameters (e.g., Quiñones et al. 2015). It is of particular interest to develop a methodology capable of extracting information of a system under changing oscillations and determine parameters that characterize the structural behavior. Usually, this is accomplished by the identification of natural frequencies of the system at different stages that can later be related to changes in stiffness (e.g., Curadelli et al. 2008;Aguirre et al. 2013;Aguirre and Montejo 2014). However, these methods are limited due to the low sensitivity of natural frequencies to localized irregularities/damage (Curadelli et al. 2008) and the saturation of the frequency shifts at large levels of ductility demand (Aguirre and Montejo 2014). In the case of damping, some experimental tests and simulations have shown that damping could be more sensible to the occurrence of damage in structural systems (Curadelli et al. 2008), while others have been found no conclusive correlation between damping ratios and the level of damage (Aguirre and Montejo 2014). Nevertheless, robust estimation of an instantaneous damping parameter could play a key role in any structural health monitoring system. Since identification of instantaneous damping relies on the estimation accuracy of both instantaneous amplitude and instantaneous frequency, the objective of this research is to propose an alternative approach to improve the extraction of these instantaneous parameters from stationary and non-stationary decaying signals. Moreover, it is shown that direct application of current methodologies based on Hilbert and/or Wavelet Transforms exhibits serious limitations due to round-up errors, edge effects, and noise contamination in the signal. To evaluate the performance of the methodologies, signals with different damping characteristics are tested for non-stationary conditions and the presence of noise contamination.

Identification via Hilbert Transform
Damping identification in a monocomponent signal x(t) can be accomplished using the Gabor analytic signal z(t) obtained from the Hilbert Transform (HT) of x(t) The analytic signal allows identification of the instant amplitude a(t) and instant phase h(t): which can be rewritten as The instant frequency (IF) is then defined as the time derivative of the phase function. To avoid discontinuities in the derivative, IF can be estimated from the phase of the derivative (Feldman 2011) as shown in Eq. 4 for a discrete signal.

Hilbert Transform limitations
Errors in calculating the discrete Hilbert Transform (HT) of a signal arise from three main reasons: (1) end effects: due to the finite duration of the signal, (2) discretization: through the sampling process, and (3) the Bedrosian identity is not satisfied because the exponential function has an unbounded spectrum. In Meissner (2012), a general error model was used to estimate the Hilbert transform of a signal in terms of the number of samples n as where the terms e 1 ½n correspond to the multiplication source of error causing erroneous oscillations of the instantaneous parameters and e 2 ½n corresponds to the addition source contributing to the non-oscillating component of the error. The error of the transform is better illustrated in the analysis of a chirp signal with increasing damping which represents a critical analysis signal with aggravated errors. The results of the application of the HT to a 20-s down-chirp signal with frequency varying from 3 to 2 Hz and linearly increasing damping are shown in Fig. 1, for the sake of clarity only the last 10 s of the signal is shown. In this case, the complementary signal (imaginary part) goes out of alignment from the real component and the accuracy of the analysis is severely affected from half the length of the signal and forth. Down-sampled versions of the imaginary part and instant amplitude are also displayed so that the oscillating behavior can be appreciated.

Improving the application of the Hilbert Transform
Even at small damping levels, the relative amplitudes in the signal affect the quality of the reconstructed analytic signal and the estimation of the instantaneous parameters. An alternative to reduce the relative amplitudes is to analyze the signal as individual segments of reduced duration. Figure 2a shows the results obtained when the HT is applied to the complete length of the signal previously presented in Fig. 1. In Fig. 2b, the HT procedure is applied only to the data in the time interval [18,20] seconds. In this case by analyzing only a portion of the signal, the results are improved significantly with the resulting components of the analytic signal aligned at a relative phase. Based on these results, the quality of the analytic signal obtained via HT can be improved by implementing a windowing procedure (Fig. 3): 1. Divide the signal x(t) into M segments of equal length N. 2. Identify additional M-1 segments of the same length N, starting from the center point of the first segment identified in step 1. The segments in this step will overlap by half their length the segments from step 1. 3. Apply the Hilbert transform operations to all of the segments from steps 1 and 2 and obtain the corresponding analytic signals. 4. The analytic signal for the whole data set is formed by putting together the results from each independent segment as shown in Fig. 3 (gray segments). Alternatively, extract first the desired instantaneous parameter from the analytic signal of each independent segment and then construct the whole time history of the parameter.
As is any windowing scheme, the length of the window will limit the range of the frequencies that can be identified and should be accounted for at the time of defining the number of segments to use. Moreover, most of the edge effects are ameliorate internally in the process of intercalation and substitution of segments. Nonetheless, when the analytic signal is generated for each segment, Gibbs effects cause a small degree of incompatibility in the intercalated segments comprising the final analytic signal. Overall, a better behaved analytic signal is obtained at the expense of a localized distribution of error introduced by the combination of segments. To ameliorate this effect, it is preferable to extract the desired instantaneous parameter from the analytic signal of each independent segment and then construct the whole time history of the parameter as shown in Fig. 3. One of the advantages of applying this method is that for signals with high frequency discontinuities, these discontinuities are isolated and its effects will appear only at the segment in which the discontinuity is located (as opposed to be distributed over the whole signal in the traditional approach). Figure 4 shows the results for the analysis of the down-chirp signal from 3 to 2 Hz with linear damping ratio from 0.5 to 5 %. Notice that resulting instantaneous parameters behave much better than the results obtained using the traditional approach (Fig. 1).

Identification via continuous wavelet transform
The continuous wavelet transform (CWT) has been proved to be an effective tool for the analysis of multi-degree of freedom systems due to its high resistance to noise-polluted signals. While the CWT capabilities for system identification in signals with relative low signal to noise ratios have been demonstrated before (e.g., Boltežar and Slavič 2004), the presence of the edge effects in the wavelet analysis represents a serious limitation in the identification of damping. In this work, we propose an alternative algorithm scheme to improve the instantaneous dynamic parameter identification capabilities of the CWT. Figure 5 describes how the CWT is usually used for instantaneous frequency and instantaneous amplitude identification. When the CWT is applied to the signal using a complex wavelet form like Complex Morlet Wavelet, the result is a matrix of complex coefficients with each coefficient providing information about the signal at a specific time and scale (which can be translated into frequency). Figure 5a shows a 2D representation (wavelet map) of the CWT coefficients using the Complex Morlet Wavelet for a linear down-chirp signal (3-2 Hz) with damping ratio increasing linearly from 0.5 to 5 %. From this map, the instantaneous frequency can be estimated by extracting the ridge (or skeleton). A complex signal proportional to the analytic signal can be constructed using the complex wavelet coefficients along the ridge (Fig. 5b) and the instantaneous amplitude to be used for damping identification can be generated (Fig. 5c). Note that edge effects cause the instantaneous amplitude envelope to displace out from the actual pattern (Fig. 5c). The results deteriorate if the signal amplitude decreases at a faster rate, for example, by increasing the damping ratio and/or the frequencies in the signal (as higher frequencies damp out faster). Figures 5d, f show the results obtained if the frequencies in the chirp previously analyzed are modified to 10-1 Hz. In this case, the analytic signal components obtained from the wavelet ridges are not properly correlated and consequently the behavior of the instantaneous parameters is erratic. The idea behind the edge effects is that the existing nonsymmetry of the wavelet function when the analysis is performed at the boundaries of the signal causes error in the phase detected and drives the envelope out of proportion. In order to ameliorate edge effects, it is common to perform a ''padding'' at the beginning and end of the signal. Proposed schemes include zero padding, value padding, decay padding, repeating the signal, reflecting the signal (e.g., Kijewski and Kareem 2002;Addison 2002). Other more unconventional techniques include no-linear extrapolation of the signal (Zheng et al. 2000) and modifications of the wavelet function when applied at the signal borders (Boltežar and Slavič 2004).
Similar to the proposed for the HT approach, an alternative to ameliorate these effects is to apply the CWT to localized segments of the signal so that the analysis is performed to portions of the signal with relatively low variations in amplitude. To avoid repeated edge effects within the segments, the proposed method consists in applying the CWT analysis to overlapping windows of constant length that are shifted a signal time step at the time. Because of the Gaussian-type amplitude envelope of the Complex Morlet wavelet, when the transform is applied the coefficient located at the center of the window is the one with the highest probability of preserving the correct instantaneous characteristics. Therefore, this is the only value kept for each window. Figure 6a, b show the analysis of the instantaneous parameters of the signals for a CWT overlapping scheme using 1 s segments. It is seen that the instantaneous frequency is properly identified; however, the instantaneous amplitude does not correlate with the actual values. Following Heisenberg uncertainty principle, there is a correlation in the localization of frequency and time. In particular, a shorter window length will result in a more accurate determination of the instantaneous amplitude parameter to the expense of decreasing the accuracy in the identification of the instant frequency. For instance, Fig. 6d and e show the analysis performed by analyzing segments of 0.5 s. In this case, the scatter in the timefrequency projection is increased while the instantaneous amplitude is extracted with more accuracy. Therefore, independent identification of the instant frequency and

Identification of instantaneous damping ratio
Typically, an average damping ratio (rather than an instantaneous damping ratio) is estimated using Eq. 6 (e.g., Montejo and Vidot-Vega 2012) where where x n is the system natural frequency, x d is the damped natural frequency, and f is the damping ratio. To obtain the average damping ratio from Eq. 6, the value of fx n can be estimated as the slope of the instant amplitude previously linearized by taking the natural logarithm. If the frequency of the system is not changing with time, x d can be estimated as the average value of the IF vector in Eq. 4 or from the Fourier spectrum of the signal. Additional methodologies for average damping estimation are presented in Boltežar et al. (2003).
If it is expected that the damping will be changing with time, for example due to changes in the physical properties of the system, an estimate of the instantaneous damping (rather than an average value over time) is desired. To obtain the instantaneous damping coefficient, an equation was proposed by Feldman (2011) that depends on the instantaneous parameters, however it relies in the first and second derivatives, which introduces discontinuities and may require post-processing. An alternative estimation is proposed here to estimate the instantaneous damping ratio where c t ð Þ ¼ Equation 8 allows estimation of the damping ratio as function of time without introducing or spreading discontinuities other than those of the instantaneous parameters. The value of the initial amplitude (IA o ) is critical and will determine the shape of the damping pattern obtained and whether it converges to the actual damping values. Figure 7 shows the envelope of the damping ratio for different values of IA o for a decaying signal with initial amplitude of 4 and a constant damping ratio of 1 %. It is seen that if the value of the initial amplitude differs from the actual value, the obtained damping values and behavior diverge from the target especially at the beginning of the signal.

Comparison between wavelets and Hilbert approaches and the effect of noise
A comparison of the results obtained with each on the two approaches proposed and its robustness to noise contamination is presented in this section. Figure 8 shows the damping results obtained using Eqs. 8 and 9 for a downchirp signal (5-1 Hz) with linear damping ratio varying linearly from 0.5 to 5 % (no noise added). The figures on the left column show the results obtained using the HT approach using 5 segments and the right column figures correspond to the CWT approach using a window length of 0.5 s. The dashed blue lines denote the target values and the continuous black lines the results obtained. It is seen that for the case of a clean signal, both approaches exhibited acceptable results. The effect of noise is examined in Fig. 9. The signal analyzed is the same analyzed in Fig. 8 but this time Gaussian noise with a SNR of 10 is added. It is seen that instantaneous parameters estimated using the CWT approach are considerably more stable and closer to the target values. For both approaches, the value of IA o to use in Eq. 9 is approximated as the maximum amplitude in the instant amplitude vector.

Conclusions
Current limitations in the identification of instantaneous damping from the free decay response of non-linear systems were identified. Improvements to the current Hilbert and continuous wavelet transform-based approaches to generate the system analytic signal were proposed so that they are suitable for identification of instantaneous parameters in vibrational signals from systems where both, damping and frequency, are changing with time. An alternative procedure to identify the instantaneous damping history from the analytic signal was also presented. It was shown that the improved approaches offer acceptable results when the signal being analyzed are clean. However, when the signal is contaminated with noise, the CWT approach exhibited better results. This can be attributed to the inherent filtering nature of the operations performed: from the resulting wavelet map, the noise content is scattered throughout the map and only the time-frequency portion (ridge) of the map that describes the system characteristics is used.