A new instantaneous frequency extraction method for nonstationary response signals in civil engineering structures

Abstract

Civil engineering structures in essence belong to time-varying and nonlinear structures and the resultant dynamic responses are usually nonstationary signals. A new method is proposed to extract instantaneous frequency of such nonstationary signals. The new method combines the extended analytical mode decomposition, normalization scheme, and improved arccosine function algorithm. In this method, a multicomponent response signal is first decomposed into several mono-components signals by the extended analytical mode decomposition theorem, and then the extracted mono-components are demodulated into purely frequency modulated signals by a repeated normalization scheme. After that, the improved arccosine function is presented to compute the instantaneous phase and instantaneous frequency. A numerical example and an experimental study of a cable with time-varying tension force are provided to illustrate the working of the method. The results demonstrate that the proposed method can extract instantaneous frequencies from multicomponent signals more accurately than Hilbert–Huang transform.

Introduction

Practical civil engineering structures essentially belong to time-varying and nonlinear structures and the resultant dynamic responses are often nonstationary due to the variation of dynamic properties under service loads or extreme loads. However, traditional signal processing techniques, including time- and frequency-domain analyses, are based on the assumption that the response signals are stationary and linear. ¹ For example, the Fourier transform can only give meaningful interpretation to linear and stationary signals, but its application to nonstationary signals of time-varying or nonlinear structures is often problematical. ² In contrast, time–frequency analysis (TFA) method is potential to process nonstationary signals and characterize their time-varying features. Various TFA methods, including short-time Fourier transform, ³ quadratic distributions, ⁴ empirical mode decomposition (EMD)-based Hilbert–Huang transform (HHT),^5,6 ensemble empirical modal decomposition (EEMD), ⁷ and continuous wavelet transform,^8,9 have been proposed in recent years to deal with nonstationary signals. However, for some civil engineering structures, the response signals may be composed of closely spaced or low-frequency components. In such cases, some classic TFA methods such as Hilbert transform (HT) cannot work well. To address this issue, it is essential to introduce effective methods for signal decomposition. ¹⁰ EMD is designed to extract intrinsic mode functions from response signals adaptively by cubic spline interpolation. However, EMD-based methods are empirical methods and still face several challenges in some engineering applications. For example, EMD is difficult to decompose closely spaced frequency components. ¹¹ By injecting a generated white noise into the original signal, EEMD may be able to separate closely spaced components but the corresponding calculation amount increases dramatically. In addition, the signal decomposition results using EEMD are sensitive to parameter selection, even pseudo-components may be produced if improper parameters are selected. ¹²

The wavelet transform as an advanced TFA technique is suitable for nonstationary signal processing and it has received wide acceptance in the dynamic signal analysis of time-varying structures in recent years.^13,14 However, Kijewski and Kareem ¹⁵ argued that civil engineering structures usually possess low-frequency motions and thus require finer frequency resolution. Moreover, parameter selection, which remained not well resolved, is critical for wavelet transform and it may significantly affect instantaneous frequency (IF) extraction results. Recently, Daubechies et al. ¹⁶ proposed synchrosqueezing wavelet transform (SWT) to reshape time frequency curves and reconstruct individual mono-components. Modified algorithms and successful applications about SWT can be found in Thakur et al., ¹⁷ Li and Liang, ¹⁸ and Cao et al. ¹⁹ It should be noted here that there is an assumption to hold by default for the SWT, that is the components of response signals should be asymptotic and well separated in frequency domain. However, some civil engineering structures such as large span cable-stayed bridges with dense stayed cables may break the assumption of well-separated frequencies and asymptotic signals. To address this issue, a new signal decomposition theorem called analytical mode decomposition (AMD) ²⁰ was presented and on this basis, the author of this paper and his cooperators extends AMD theorem to solve the problem of modal mixing. ²¹ Numerical simulation results indicate that each frequency-modulated (FM) individual component between any two bisecting frequencies can be analytically extracted by extended AMD.

However, the mono-components decomposed by the extended AMD are generally amplitude-modulated (AM) and frequency-modulated (FM) functions, and the amplitude change rate of the mono-component signal is usually not much slower than its phase change rate. Therefore, the significance of the demodulation process is highlighted. The local mean decomposition (LMD) method developed by Smith ²² is a self-adaptive algorithm to decompose AM–FM signals into a small set of product functions, each of which is the product of an envelope signal and a FM signal. The LMD avoids the limitation of negative frequencies occurring in HT, but end effect and mode mixing problems remain two main limitations in LMD, which can be alleviated by proper parameter selection with regard to boundary condition, envelope estimation, and sifting stopping criterion. ²³ Huang et al. ¹⁰ proposed a normalization scheme to demodulate signals through empirical amplitude method. The normalization scheme is exactly the process of demodulation based on iterative applications of cubic spline fitting and thus the AM and FM parts of the decomposed components are separated uniquely. With this repeated demodulation process, we can conduct HT or other TFA methods on the FM part alone, avoiding the difficulty stated in the Bedrosian theorem and the Nuttall theorems. ² The normalization scheme not only successfully extract FM parts of original signals but also converge very fast in the repeated demodulation process. ¹⁰ Nevertheless, the normalization process could cause some deformation of the original signal, but the amount of the deformation is negligible, for there are rigid controlling points for the periodicity provided by the zero-crossing points in addition to the extrema.

The next step is to extract IFs from FM signals. The most common methods for calculating instantaneous phase or IF include HT, Teager energy operator (TEO), zero-crossing, and arccosine function method.^2,10 Among them, HT is the most popular and direct method to extract IF. However, the IF extracted by HT may be negative at times, which is extremely hard to interpret. In addition, HT is very sensitive to random noises and hence leads to an IF curve with blurred and distorted lines. The TEO is a method to compute IF totally based on differentiations. A distinct advantage of the TEO is its excellent localization property and no integral transform needed as in HT. The shortcoming of this method is that TEO is based on a linear model for a single harmonic component only. Therefore, the approximation produced by the TEO will deteriorate and even break down when either the amplitude is a function of time or the waveforms have any interweave modulations or harmonic distortions. ¹⁰ The zero-crossing method has long been used to compute the mean frequency for narrow band signals. But unfortunately, the calculation results are relatively crude for the frequency should be constant over the period between two zero-crossings. The arccosine function is potential to calculate the phase from FM signals, and then the IF can be expressed as the derivative of the instantaneous phase with respect to time. Although the arccosine function is able to calculate the phase successfully, its range is strictly restricted to $\left[0,\text{π}\right]$ . Therefore, the arccosine function cannot be directly applied to the purely FM signal whose phase range actually is the set of real numbers. ²⁴ To address this issue, a new method called improved arccosine function is proposed in this paper to solve the instantaneous phase and IF more accurately based on the division of intervals along time axis. The main limitation of improved arccosine function method is the distortion of the instantaneous phase, which can be greatly reduced by a moving average algorithm. ²³ With the moving average method, the instantaneous phase curves beforehand extracted by the improved arccosine function are smoothed without great distortions.

The contribution of this paper is to propose a new IF extraction method, which is different from the traditional HHT. The proposed new method, which combines the extended AMD, normalization scheme, and improved arccosine function, is composed of three steps accordingly. First, the extended AMD is employed to decompose a response signal, which may include closely spaced or low-frequency components, into several mono-components. Second, the normalization scheme is used to demodulate the decomposed mono-component signal and thus the AM and FM parts are completely separated. Finally, the improved arccosine function is presented to compute the instantaneous phase. To reduce the distortion of the instantaneous phase caused by segmentation error, a moving average algorithm is introduced to smooth the target instantaneous phase curves. Then, the IF can be solved by taking a derivative of instantaneous phase with respect to time. The main difference between the proposed method and HHT is that the decomposition is not empirical but analytical. In addition, a major highlight of the new method is the presentation of the improved arccosine function.

Theory and methodology

The method proposed in this paper combines the techniques of extended AMD, normalization scheme, and improved arccosine function. It allows to decompose and demodulate multicomponent signals composed of closely spaced frequency or low-frequency components into mono-components with separated AM and FM parts, whereby the improved arccosine function can be used to achieve good time–frequency resolutions. The flowchart of this new method is illustrated in Figure 1.

Figure 1.

The flowchart of the new method.

The extended AMD theorem is first introduced to decompose a multicomponent response signal into several mono-components by selecting proper bisecting frequencies. A normalization scheme is then employed to demodulate the decomposed mono-component signals into purely FM. After that, the improved arccosine function is presented to characterize the IFs of the purely FM signals with a better time–frequency resolution.

Signal decomposition by the extended AMD theorem

For nonstationary AM–FM signals, the frequency of each individual component varies with time. Instead of selecting constant bisecting frequencies, time-varying bisecting frequencies estimated from a preliminary SWT are selected, and each FM component between any two bisecting frequencies can then be analytically extracted by the extended AMD theorem.

Let x(t) denote a real response signal of r individual components with frequencies: ${\omega }_1(t),{\omega }_2(t),\dots ,{\omega }_r(t)$ , all positive, in L₂(−∞, +∞) of the real time variable t. It can be decomposed into r signals $x_i^{\left(d\right)}(t)$ (i = 1, 2, … , r) whose frequency ranges satisfy $|{\omega }_1\left(t\right)|<{\omega }_{b1}(t)$ , ${\omega }_{b1}(t)<|{\omega }_2\left(t\right)|<{\omega }_{b2}(t),\dots$ , ${\omega }_{b\left(r-2\right)}(t)<|{\omega }_{r-1}\left(t\right)|<{\omega }_{b\left(r-1\right)}(t)$ , and ${\omega }_{b\left(r-1\right)}(t)<|{\omega }_r(t)|$ . That is

$$x\left(t\right)={\displaystyle \sum }_{i=1}^r{x}_i^{(d)}(t)$$

(1)

Here, ${\omega }_i(t)$ is a frequency variable for the decomposed signal $x_i^{(d)}(t)$ , and ${\omega }_{bi}(t)\in \left({\omega }_i(t),{\omega }_{i+1}(t)\right)(i=1,2,\dots ,r-1)$ is time-varying bisecting frequencies. Then, each individual signal can be determined by

$$x_1^{(d)}=s_1\left(t\right),\dots ,x_i^{(d)}\left(t\right)=s_i\left(t\right)-s_{i-1}\left(t\right),\dots ,x_r^{\left(d\right)}(t)=x\left(t\right)-s_{r-1}(t)$$

(2)

and

$$s_i\left(t\right)=\sin [{\theta }_{bi}(t)]\text{H}\{x\left(t\right)\text{cos}[{\theta }_{bi}\left(t\right)]\}-\cos [{\theta }_{bi}\left(\text{t}\right)]\text{H}\{x\left(t\right)\sin \left[{\theta }_{bi}\left(t\right)\right]\}\hspace{1em}\left(i=1,2,\dots ,r-1\right)$$

(3)

where $\text{H}\left\{\right\}$ returns the HT the function inside the square bracket and ${\theta }_{bi}\left(t\right)={\int }_{-\infty }^t{\omega }_{bi}(\tau )\text{d}\tau$ is the phase function of the $i$ th bisecting frequency. ${\omega }_i(t)$ and ${\omega }_{bi}(t)$ can be estimated from a previous SWT $T_x\left({\omega }_l,b\right)$ .^16,17 More details about the proofs and numerical studies of the extended AMD theorem can be found in Wang et al. ²¹

Demodulation by a normalization scheme

For a mono-component signal $x_1\left(t\right)$ decomposed by the extended AMD theorem, the analytic signal of $x_1\left(t\right)$ can be defined as a complex-valued function which has the original signal as its real part and the HT as its imaginary part

$$z\left(t\right)=x_1\left(t\right)+\text{i}\cdot \text{H}\left[x_1\left(t\right)\right]=A_1\left(t\right)e^{-\text{i}{\varphi }_1\left(t\right)}$$

(4)

where $\text{H}[\cdot ]$ represents the HT of a signal in the square bracket, $\text{i}=\sqrt{-1}$ , $A_1\left(t\right)=\sqrt{{\left[x_1\left(t\right)\right]}^2+{\left\{\text{H}\left[x_1\left(t\right)\right]\right\}}^2}$ and ${\varphi }_1\left(t\right)=\text{arctan}\left[\text{H}\left[x_1\left(t\right)\right]/x_1\left(t\right)\right]$ are instantaneous amplitude and instantaneous phase, respectively. Then, $x_1\left(t\right)$ can be expressed in equation (5) through the definition of the analytic signal

$$x_1\left(t\right)=A_1\left(t\right)\text{cos}{\varphi }_1\left(t\right)$$

(5)

Supposing the Bedrosian theorem is satisfied, the HT for the right side of equation (5) can be rewritten as

$$\text{H}\left[A_1\left(t\right)\text{cos}{\varphi }_1\left(t\right)\right]=A_1\left(t\right)\text{H}\left[\text{cos}{\varphi }_1\left(t\right)\right]$$

(6)

Actually equation (6) can hold if the amplitude varies so slowly that the frequency spectra of the envelope and the carrier waves are disjoint according to the Bedrosian theorem. ² In other words, the signal to be analyzed should be asymptotic or purely FM, but this assumption is usually broken for response signals in practical civil engineering structures.

To address this issue, the signal normalization scheme is used to demodulate the extracted component signals into AM and FM parts uniquely. First, we calculate absolute values of the decomposed mono-component signal and then identify all extrema of these absolute values, followed by defining the envelope with a spline through all these extrema and designating it as $e_1\left(t\right)$ . As such, the normalization $d_1\left(t\right)$ can be expressed as

$$d_1\left(t\right)=x_1\left(t\right)/e_1\left(t\right)$$

(7)

However, $e_1\left(t\right)$ could have values smaller than the decomposed signal $x_1\left(t\right)$ at some temporal locations and cause the normalized function to have amplitude greater than unity when the fitted envelope $e_1\left(t\right)$ based on the extrema has cuts through $x_1\left(t\right)$ . These conditions may be rare but are not impossible. Therefore, the normalization should be further conducted recursively as equation (8)

$$d_2\left(t\right)=d_1\left(t\right)/e_2\left(t\right),\dots ,d_m\left(t\right)=d_{m-1}\left(t\right)/e_m\left(t\right)$$

(8)

This iteration process does not stop until the normalized function values are all unity. Set m as the number of total iterations, and the final demodulated FM part $D_1\left(t\right)$ and AM part $A_1\left(t\right)$ of $x_1\left(t\right)$ are given by

$$\{\begin{array}{c}{D}_1\left(t\right)=d_m\left(t\right)\\ A_1\left(t\right)=x_1\left(t\right)/D_1\left(t\right)\end{array}$$

(9)

So far the demodulation process is successfully realized and the FM signal separated from mono-component signal $x_1\left(t\right)$ can be used in the next step.

IF extraction by improved arccosine function

One way to express the nonstationarity is to find instantaneous phase and IF using HT. With HT, the instantaneous phase is solved by the definition of analytical signal expressed in equation (4). Then, the IF can be calculated by taking the derivative of instantaneous phase with respect to time. However, a purely oscillatory function with a zero reference level is a necessary condition for the abovementioned HT method. Otherwise, the distortion of phase function will happen and thus cause the deviation of the extracted IF from its true value. ⁵ Moreover, the sensitivity to noises limits its practical application. An alternative way to compute the phase from the FM signal is the arccosine function. The arccosine function can calculate the phase successfully, but its range is strictly restricted to $\left[0,\text{π}\right]$ . Therefore, the arccosine function cannot be applied to the purely FM signal whose phase range actually is a set of real numbers. As such, we present a new method called improved arccosine function based on the division of intervals along time axis to solve the instantaneous phase more accurately than the standard arccosine function.

Typically the time domain graph of a purely FM signal $x_1\left(t\right)$ is composed of several waves as shown in Figure 2, with each having a wave crest and trough. We can find two zero-crossing points in each wave by using the operator $\text{sgn}()$ .

Figure 2.

The multicomponent signal $\mathbf{x}\left(\mathbf{t}\right)=\{\begin{array}{cc}\mathbf{cos}\mathbf{2}\mathbf{\pi }\mathbf{t}& \mathbf{0}\le \mathbf{t}\le \mathbf{5}\\ \mathbf{cos}\mathbf{4}\mathbf{\pi }\mathbf{t}& \mathbf{5}<\mathbf{t}\le \mathbf{10}\end{array}$ with its zero-crossing points and extrema.

Set n as the total time points of the purely FM signal $x_1\left(t\right)$ , and define a new time series $s\left(t_i\right)$ in equation (10)

$$s\left(t_i\right)=x_1\left(t_{i-1}\right).x_1\left(t_i\right)\hspace{1em}(i=2,\dots ,n)$$

(10)

Then, the operator $\text{sgn}()$ is expressed as equation (11)

$$\text{sgn}\left(s\left(t_i\right)\right)=\{\begin{array}{cc}1,& s\left(t_i\right)>0\\ 0,& s\left(t_i\right)=0\\ -1,& s\left(t_i\right)<0\end{array}$$

(11)

If $\text{sgn}\left(s\left(t_i\right)\right)=-1$ , it implies that there is at least one time instance between two points $t_i$ and $t_{i-1}$ at which the value of $x_1\left(t\right)$ equals zero. Then the time point $t_i$ can be approximated as zero-crossing point with i as the index, provided that the sampling frequency for the signal $x_1\left(t\right)$ is large enough. In other words, the point $t_i$ can be directly defined as the zero-crossing point when $\text{sgn}\left(s\left(t_i\right)\right)=0$ . In this way, we can find all p zero-crossing points and define the indexes of them as $z_1,z_2,\dots ,z_p$ . Then the time domain graph of a purely FM signal $x_1\left(t\right)$ can be divided into p + 1 intervals named as [ $0,t_{z_1}$ ], [ $t_{z_1},t_{z_2}$ ], …, [ $t_{z_p},t_n$ ], and each individual maximum/minimum extremum in the intervals can be extracted, with the indexes of the extrema are named as $c_1$ , $c_2$ , …, $c_p$ , $c_{p+1}$ . By using these extracted extrema, the time axis is separated into p + 1 intervals and defined as [ $0,t_{c_1}$ ], [ $t_{c_1},t_{c_2}$ ], … , [ $t_{c_{p-1}},t_{c_p}$ ], [ $t_{c_p},t_{c_{p+1}}$ ].

Taking the signal $x_1\left(t\right)=\{\begin{array}{cc}\text{cos}2\text{π}t& 0\le t\le 5\\ \text{cos}4\text{π}t& 5<t\le 10\end{array}$ for example, we can find two extrema points and one zero-crossing point at any time interval [ $t_{c_i},t_{c_{i+1}}$ ] after assigning the operator $\text{sgn}(),$ which is shown in Figure 2. Assuming $x_1\left(t\right)$ is monotonic in each individual interval [ $t_{c_i},t_{c_{i+1}}$ ], we can draw a conclusion that the phases of $x_1\left(t\right)$ in the first interval fall in the range of $\left[0,\text{ π}\right]$ . Likewise, the phases of $x_1\left(t\right)$ in the second interval are located in $\left[\text{π},2\text{π}\right]$ . By analogy, the phases of $x_1\left(t\right)$ in the ith interval will definitely be $\left[(i-1)\text{π},i\text{π}\right]$ . Therefore, the instantaneous phase ${\varphi }_1\left(t\right)$ of $x_1\left(t\right)$ can be accurately computed as follows

$${\varphi }_1\left(t\right)=\{\begin{array}{cc}\text{arccos}\left[x_1\left(t\right)\right]& t\in \left[0,t_{z_1}\right]\\ {(-1)}^{i-1}\cdot \text{arccos}\left[x_1\left(t\right)\right]+2\text{π}\cdot \text{floor}\left(\frac{i}{2}\right)& t\in \left[t_{z_i},t_{z_{i+1}}\right]\end{array}$$

(12)

where $\text{floor}(\cdot )$ represents the element inside the parenthesis to the nearest integer toward minus infinity.

However, the distortions of the phases usually occur at the locations near extrema points. The reason for this phenomenon is that the extrema points defined in the improved arccosine function algorithm can only be infinitely close to the true extrema points due to the discretization of signals and thus the occurrence of distortions near the center of the extrema is not a surprise. In addition, noises have an impact on the calculation of instantaneous phase, but the effect is small because most high-frequency noises are already removed by the extended AMD filter. In order to solve the problem of phase distortions, it is necessary to find a useful mathematical tool to smooth the whole instantaneous phase curve. Moving average method is one of useful methods for smoothing the computed instantaneous phases. A simple moving average method is actually a common average of n consecutive data points in time series with each point equally weighted. In this paper, only five-point moving average method is used and the formula is defined as follows

$$f_a\left(t_i\right)=\frac{f\left(t_{i-2}\right)+f\left(t_{i-1}\right)+f\left(t_i\right)+f\left(t_{i+1}\right)+f\left(t_{i+2}\right)}{5}$$

(13)

where $f\left(t_i\right)$ stands for the value of the frequency at $t_i$ and $f_a\left(t_i\right)$ represents the averaged or smoothed frequency at $t_i$ . It is noted that the smoothed frequencies at the first two and the last two points cannot be directly calculated by equation (13). One possibility is to update the old value at abovementioned locations with the counterparts at neighboring locations. By using the moving average method, the computation stability and identification accuracy of instantaneous phases along the whole time axis can be greatly enhanced.

Thus, the IF can be further computed by

$$f\left(t\right)=\frac{\omega \left(t\right)}{2\text{π}}=\frac{1}{2\text{π}}\frac{\text{d}{\varphi }_1\left(t\right)}{\text{d}t}$$

(14)

in which $\omega \left(t\right)$ and ${\varphi }_1\left(t\right)$ represent the circular frequency and the instantaneous phase of the decomposed component signal, respectively. Since the variables $f\left(t\right)$ and ${\varphi }_1\left(t\right)$ are often discretized in practical engineering, the discrete version of equation (14) is expressed as

$$f\left(t_i\right)=\frac{1}{2\text{π}}\frac{{\varphi }_1\left(t_i\right)-{\varphi }_1\left(t_{i-1}\right)}{t_i-t_{i-1}}$$

(15)

Numerical examples

In this section, a numerical simulation case is employed to verify the effectiveness and accuracy of the proposed method. In addition, an index of accuracy (IA) is presented to quantify the accuracy of the IF identification.

A multicomponent signal with two closely spaced AM–FM components is considered

$$y\left(t\right)=y_1\left(t\right)+y_2(t)$$

(16)

where $y_1\left(t\right)=[3+\cos \left(\text{π}t\right)]\text{cos}[8\text{π}t+2\sin \left(\text{π}t\right)]$ and $y_2\left(\text{t}\right)=[4+\cos \left(\text{πt}\right)]\text{cos}[13\text{πt}+\sin \left(2\text{πt}\right)]$ . The theoretical frequencies of the two components are closely spaced and they are expressed as $f_1=4+\text{cos}\left(\text{πt}\right)$ Hz and $f_2=6.5+\text{cos}(\text{πt})$ Hz, respectively. The sampling frequency is 100 Hz and the sampling time is 10 s. A simulated 20% Gaussian white noise is added to consider the influence of noises. The noise intensity is defined by the following signal-to-noise ratio (SNR)

$$\text{SNR}=10{\text{log}}_{10}\frac{A_{\text{signal}}^2}{A_{\text{noise}}^2}\left(\text{in dB}\right)$$

(17)

where $A_{\text{signal}}$ and $A_{\text{noise}}$ refer to the root mean square value of the signal and the noise, respectively, and the noise level means the ratio of $A_{\text{noise}}^2$ to $A_{\text{signal}}^2$ .

The waveform of the noisy multicomponent signal to be analyzed is plotted in Figure 3. The SWT with Morlet transform as parent function is first performed on the noisy signal to obtain the wavelet scalogram of $y(t)$ as presented in Figure 4. By selecting 5 Hz as the bisecting frequency, two individual components can be extracted by the extended AMD, which are shown in Figure 5. Then the normalization scheme is used to separate the decomposed mono-components into AM and FM parts until the maximum is close to 1. The demodulation results are also illustrated in Figure 5. One can see from Figure 5 that the decomposed mono-components extracted by the extended AMD agree with theoretical results well while the counterparts extracted by EMD do not. Moreover, the periodic change of waveform of the purely FM signals obtained by the normalization scheme is also in good accordance with that of theoretical counterparts. Finally, the instantaneous phases and IFs are computed by the improved arccosine function and then smoothed by the moving average method. Figure 6 presents the results of IF extraction using the new method and HHT, which indicates that the new method can provide a more accurate IF curve than HHT.

Figure 3.

The simulated AM–FM multicomponent signal.

Figure 4.

Wavelet scalogram of the AM–FM multicomponent signal.

Figure 5.

Extracted mono-components: (a) y₁ and (b) y₂. AMD: analytical mode decomposition; EMD: empirical mode decomposition.

Figure 6.

IFs extracted by the new method and HHT: (a) y₁ and (b) y₂. HHT: Hilbert–Huang transform.

To quantify the accuracy of the IF identification, an IA is defined as the root-mean-squared value over the total time duration of the IF and expressed as equation (18)

$$\text{IA}=\sqrt{\frac{{{\displaystyle \int }}_0^T{\left[f_d(t)-f_e(t)\right]}^2}{{{\displaystyle \int }}_0^T{\left[f_e(t)\right]}^2}}$$

(18)

where $f_d(t)$ is the extracted IF, and $f_e(t)$ is the exact/theoretical IF. ${\text{IA}}_1$ and ${\text{IA}}_2$ represent the IF identification accuracies of components $y_1(t)$ and $y_2(t)$ , respectively, The smaller the value of $\text{IA}$ , the better accuracy of IF extraction.

According to equation (18), the values of IA using HHT and the new method are computed and then presented in Table 1. We can see from Table 1 that the value of IA using the new method is much smaller than that of the HHT, which verifies the superiority of the new method, especially when response signals are contaminated by random noises.

Table 1.

IA of IF extraction of the noisy signal with AM–FM components (%).

Methods	New method	HHT
${\text{IA}}_1$	5.1	113.4
${\text{IA}}_2$	6.7	455.7

AM: amplitude modulated; FM: frequency modulated; HHT: Hilbert–Huang transform; IA: index of accuracy; IF: instantaneous frequency.

Experimental case study

IF identification of a cable with linearly varying tension forces

To further validate the accuracy of the proposed method, a cable composed of pieces of 7Φ5 steel wires is considered. The cable has elastic modulus E = 1.95 × 10⁵ MPa, area of cross section A = 1.374 × 10⁴ m², and density of unit length q = 1.1 kg/m. The stiffness of the cable is changed due to the applied time-varying tension force so that the natural frequency of the cable is time dependent. The cable is fixed at one end, and the other end is connected to an MTS loading system. The total length of cable between two ends is 4.55 m. A accelerometer is installed vertically at the midpoint of the cable. The test setup is shown in Figure 7.

Figure 7.

The cable test setup.

During the cable test, an initial constant pretension force was first applied by MTS load actuator to the cable. When the test setup and data collection were ready, the cable tension force was changed continuously using the MTS load system. At the same time, the impact hammer was used to generate free vibration and the vertical acceleration responses were recorded at a sampling frequency of 600 Hz. The tension force with linear or sinusoidal change was considered during the test. The theoretical IFs of the cable were obtained by solving eigenvalues and eigenvectors of vibration equations, assuming that the parameters of the cable keep invariant over a relatively short time interval, which is designated as the time frozen method. ⁹ For simplicity, only the fundamental frequency was considered at this case. The theoretical fundamental frequencies of the cable at different fixed tension forces are listed in Table 2. All these fundamental frequencies were used as theoretical IFs at different constant tension forces to compare with the extracted IFs under time-varying tension forces.

Table 2.

The theoretical fundamental frequency of the cable under different constant tension forces.

Tension (kN)	13.0	15.0	17.0	19.0	20.0	20.3	20.6	20.9	21.2
Frequency (Hz)	12.26	13.08	13.87	14.65	15.03	15.13	15.24	15.33	15.45
Tension (kN)	21.5	21.8	22.0	22.5	23.0	23.5	24.0	24.5	25.0
Frequency (Hz)	15.56	15.67	15.70	15.91	16.09	16.24	16.41	16.54	16.72
Tension (kN)	25.5	26.0	26.5	27.0	27.5	28.0	29.0
Frequency (Hz)	16.88	17.03	17.19	17.33	17.47	17.63	17.90

For the first case, the initial tension force of the cable was set to 20 kN, and then the tension force increased linearly at the rate of 1.67 kN/s using the MTS load system. The duration of data acquisition was set to be 6 s. The measured tension force and acceleration responses were shown in Figures 8 and 9, respectively.

Figure 8.

Measured cable tension forces with linear variation.

Figure 9.

Measured cable acceleration responses with linearly varying tension forces.

By performing SWT on acceleration responses with Morlet wavelet as parent function, the wavelet scalogram was obtained and displayed in Figure 10, which indicates that there were two IF trajectories in the time–frequency plane. Extended AMD was used to extract the fundamental frequency component and then the normalization scheme was introduced to demodulate it. The results of component decomposition and demodulation are presented in Figure 11, which shows that the periodic change of waveform of the FM signal obtained by the normalization scheme is in good accordance with that of theoretical counterparts. Finally, the improved arccosine function and moving average algorithm were implemented to identify the instantaneous phases and IFs of the cable and the theoretical fundamental frequencies solved by the time frozen method are shown in Table 3. Figure 12(a) presented the theoretical fundamental frequencies and the extracted IFs using the proposed new method and HHT. It can be seen from Figure 12 that the new method can extract IFs of the cable subjected to linearly increased tension forces more accurately than HHT.

Figure 10.

Wavelet scalogram of cable acceleration responses with linearly varying tension forces.

Figure 11.

Extracted mono-components of cable acceleration responses with linearly varying tension forces. AMD: analytical mode decomposition; EMD: empirical mode decomposition.

Table 3.

The theoretical fundamental frequency of the cable with linearly varying tension forces.

Time (s)	Tension (kN)	Theoretical frequency (Hz)	Time (s)	Tension (kN)	Theoretical frequency (Hz)
0	20.0	15.03	2.76	24.5	16.54
0.06	20.3	15.13	3.06	25.0	16.72
0.43	20.6	15.24	3.35	25.5	16.88
0.62	20.9	15.33	3.65	26.0	17.03
0.78	21.2	15.45	3.96	26.5	17.19
0.97	21.5	15.56	4.26	27.0	17.33
1.15	21.8	15.67	4.56	27.5	17.47
1.26	22.0	15.70	4.86	28.0	17.63
1.55	22.5	15.91	5.47	29.0	17.90
1.85	23.0	16.09	6.09	30.0	18.20
2.15	23.5	16.24	6.70	31.0	18.47
2.46	24.0	16.41

Figure 12.

Extracted IF of the cable with linearly varying cable tension forces. (a) Compared with HHT and theoretical IFs and (b) compared with theoretical IFs alone. HHT: Hilbert–Huang transform.

To compare the performance of the new method with that of HHT, the values of IA are calculated according to equation (18) and hence presented in Table 4. It can be concluded from Table 4 that the proposed new method is much better than HHT as for the IF extraction of a cable with linearly varying tension forces.

Table 4.

IA of IF extraction with linearly varying cable tension forces (%).

Methods	New method	HHT
IA	0.6	126.9

HHT: Hilbert–Huang transform; IA: index of accuracy; IF: instantaneous frequency.

IF identification of a cable with sinusoidally varying tension forces

In the second case, the initial tension force of the cable was set to 22 kN, and then the tension force varied sinusoidally, which is shown in Figure 13. The acceleration response signal with a duration of 6 s was measured at the middle of the cable and shown in Figure 14.

Figure 13.

Measured cable tension forces with sinusoidal variation.

Figure 14.

Measured cable acceleration responses with sinusoidally varying tension forces.

The SWT with Morlet wavelet as parent function was first conducted on the acceleration response signal and the corresponding wavelet scalogram is plotted in Figure 15. As shown in Figure 15(a), the measured acceleration response signal contains several components. Figure 15(b) illustrates that the IF of the fundamental frequency component oscillates sinusoidally. Again, the extended AMD was used to decompose the original multicomponent response signals into several mono-components and the normalization scheme was then used to demodulate the fundamental frequency component. The results of extended AMD and normalization scheme are displayed in Figure 16. It can be seen from Figure 16 that there is no obvious change of the periodicity of waveform before and after the normalization scheme. Finally, the improved arccosine function was employed to extract the IFs of the cable, which were then smoothed by the moving average algorithm. The theoretical fundamental modal frequencies were solved using the time frozen method and shown in Table 5. The comparison of theoretical fundamental frequencies and identification results presented in Figure 17 indicates that the identified IFs using the new method are in reasonable accordance with theoretical results. However, the classical HHT method cannot effectively extract IFs because of its inherent defects and the noise contamination. Moreover, the IA presented in Table 6 shows again that the new method provides better time–frequency resolutions than HHT on IF extraction of a cable with sinusoidal tension forces.

Figure 15.

Wavelet scalogram of cable acceleration responses with sinusoidally varying tension forces.

Figure 16.

Extracted mono-components of cable acceleration responses with sinusoidally varying tension forces. AMD: analytical mode decomposition; EMD: empirical mode decomposition.

Table 5.

The theoretical fundamental modal frequency of the cable with sinusoidally varying tension forces.

Time (s)	Tension (kN)	Theoretical frequency (Hz)	Time (s)	Tension (kN)	Theoretical frequency (Hz)
0	22	15.70	3.4	17.96	14.24
0.2	21.84	15.68	3.8	18.03	14.27
0.4	21.55	15.58	4.2	18.37	14.40
0.6	21.17	15.44	4.7	19.03	14.66
0.9	20.59	15.24	5.3	20.06	15.05
1.2	20.09	15.06	5.4	20.25	15.11
1.8	19.1	14.69	5.6	20.63	15.25
2.1	18.71	14.54	5.9	21.22	15.46
2.5	18.3	14.38	6.0	21.42	15.53
3.1	17.98	14.25

Figure 17.

Extracted IF of the cable with sinusoidally varying tension forces. (a) Compared with HHT and theoretical IFs and (b) compared with theoretical IFs alone. HHT: Hilbert–Huang transform.

Table 6.

IA of IF extraction of the cable with sinusoidal varying tension forces (%).

Methods	New method	HHT
IA	2.9	136.1

HHT: Hilbert–Huang transform; IA: index of accuracy; IF: instantaneous frequency.

Conclusions

This paper presents a new IF extraction method for nonstationary response signals in civil engineering structures. The proposed method is a combination of the extended AMD, normalization scheme, and improved arccosine function. The extended AMD is able to decompose a multicomponent signal, which may consist of closely spaced modal or low-frequency components, into mono-component signals. The normalization scheme demodulates the extracted mono-component signals into purely FM signals. Then the improved arccosine function with the moving average algorithm as a smoother is presented to enhance the time–frequency resolutions of IF curves. A numerical example of a signal with closely spaced AM–FM components and an experimental case study on a cable with time-varying tension forces are employed to verify the effectiveness and accuracy of the proposed method. The results demonstrate that the proposed new method can generate more accurate IF identification results than HHT. As the new method only needs the output response data, therefore, it has the potential to trace the instantaneous properties of time-varying structures under ambient excitations.