A comparative study of signal processing methods for structural health monitoring

3.1. Fourier transform

Fourier transform (FT) is the most widely-used method to represent the frequency domain of a time series signal. Assuming that the signal $\mathbf{x}\left(t\right)$ has finite energy, that is:

the discrete Fourier transform (DFT) is defined as:

where $N$ is the length of signal (number of sampling points), $f=k/N$ is the frequency in cycles per sample, $\omega =2\pi k/N$ is the circular frequency in radians/sample, and $i$ is the sample index (e.g. time increment). The Fourier transform of the synthetic signal (Fig. 1) is shown in Fig. 2. The estimated frequencies (the first four peaks noted by triangles in Fig. 2) are 0.33, 0.99, 1.49, and 1.74 Hz. Table 1 presents the estimated values of frequency and damping exponents for the methods studied in this research. In the results of DFT method, the first two damped frequencies are mixed as one peak at around 0.33 Hz. This is a drawback of the FT method that cannot extract all the frequencies accurately when the signal is a) noisy, b) has closely-spaced frequencies, or c) is non-stationary. Moreover, the 4th peak with a value of 1.74 Hz, is a spurious solution, another drawback of the FT. Furthermore, FT provides no information about the damping values.

3.2. Periodogram estimate of power spectral density (PSD)

The power spectral density (PSD), a measure suitable for analysis of stationary signals, is defined as the Fourier transform of the autocorrelation sequence of the sample signal. It is often estimated by a periodogram with a rectangular window, defined by magnitude squared of DFT coefficients as follows:

where $T_s$ is the sampling time interval and $f$ is frequency in Hz, $0<f\le 1/2T_s$ due to Nyquist Sampling theorem. The PSD of the sample signal is shown in Fig. 3. Table 1 presents the estimated values of the four frequencies. Similar to DFT, it represents two frequencies accurately but the first two closely-spaced modes are mixed together and indistinguishable. The 4th frequency at 1.73 is a spurious solution. As such, PSD also suffers from both problems of mixing closely-spaced frequencies and generating spurious solutions. Further, PSD shows noisy characteristics with many spurious peaks.

Fig. 3Power spectral density (PSD) of the signal in Fig. 1 with the peak frequency values of 0.31, 1.02, 1.57, and 1.73

3.3. Wavelet transform

In the past two decades wavelet transform has been used extensively as a time-frequency signal processing method in a variety of applications such as seismology [33]; structural and mechanical vibrations [34], vibration control [35], system identification [36], power engineering [37], computational neuroscience [38, 39], traffic engineering [40-43] and SHM [44]. In contrast to FT where sinusoidal functions are used, in continuous wavelet transform a base function with localized energy is used to approximate the signal. The signal is decomposed into different variations of the base function and presented as an integral of scaled and dilated wavelets. The coefficients are found from the following transformation represented by a two-dimensional matrix $\mathbf{W}$:

where $a$ and $b$ are the scale and the dilation parameters, and $\mathrm{\Psi }$ is the base wavelet function, called mother wavelet.

In discrete wavelet transform (DWT), two quadrature mirror filters, one low pass and one high pass are used in the first level of transformation. The results from high pass and low pass filters are called detail and approximation coefficients, respectively. Since each output contains half of the frequency range, they are subsampled by 2. As such, the time resolution will be only half of the original signal. This process is repeated to calculate all coefficients representing higher frequency resolutions. To improve the accuracy, wavelet packet transform (WPT) [45, 46] has been proposed where both approximation and detail signals are decomposed. Fig. 4 shows the values of the calculated WPT coefficients of the synthetic signal shown in Fig. 1 in time and frequency domains using the widely-used Daubechies wavelet of the 1st order. This figure is known as a scalogram. The estimated values of the four frequencies presented in Table 1 indicate a concentration of the frequency content at about 0.5, 1, and 1.5 Hz noted by dark rectangles and the dissipation in time, as expected, noted by lightening of the rectangles. However, the time-resolution is coarse and the closely-spaced frequencies at 0.3 and 0.5 Hz are mixed together as one and indistinguishable. The result may improve by choosing a higher frequency resolution but at the cost of a lower time resolution. Also, the wavelet type selected has some bearing on the accuracy which makes it problem-dependent or a subjective decision making process. WT and WPT are indeed effective for detecting sudden changes in a signal but less effective for extracting the contents of smoothly changing or exponentially decaying signals.

Fig. 4Scalogram of the synthetic signal shown in Fig. 1 using the wavelet packet decomposition

3.4. Empirical mode decomposition (EMD) with Hilbert spectral analysis (Hilbert-Huang transform)

The empirical mode decomposition (EMD) method is a time domain signal processing approach that decomposes a signal into Intrinsic Mode Functions (IMFs). It is used as the first step of Hilbert-Huang transform (HHT) where Hilbert spectral analysis is performed on the extracted IMFs to obtain the instantaneous frequency contents of the data [47]. HHT was proposed mainly for analysis of non-stationary nonlinear signals. Any given time domain signal $\mathbf{x}\left(t\right)$ is decomposed into $n-1$ oscillating and one non-oscillating IMFs:

The method yields all modes with no prior information about the number of signal components ($n$) which makes it attractive for non-parametric signal processing. First, all local maxima and minima of the original signal are detected. Then, two envelopes including all maxima and minima are calculated as a function of time, $\overline{\mathbf{x}}\left(t\right)$ and $\underset{\_}{\mathbf{x}}\left(t\right)$, respectively. Then, the mean of these two time functions is deducted from the original signal:

This process is repeated for the residual ${\mathbf{x}}^1\left(t\right)$ signal until the mean vector approaches a null vector within a given tolerance. The last residual sequence is the first extracted IMF, $IM{F}_1\left(t\right)$, of the signal which usually includes the highest frequency content. This process is repeated for $\mathbf{x}\left(t\right)-\sum _{i}IM{F}_i\left(t\right)$ until the last remainder has only one maximum or minimum.

Next, an imaginary signal introduced by [48] as the analytic signal is created as follows:

where $\mathbf{y}\left(t\right)$ is the Hilbert transform of the signal, $\mathbf{x}\left(t\right)$:

Hilbert transform of a signal is in fact the convolution integral of the signal and $1/\pi t$. The analytic signal as a complex signal can be expressed in terms of the magnitude and the phase angle, $\theta \left(t\right)$, where the instantaneous frequency is the time derivative of the phase angle:

Fig. 5(a) and Fig. 5(b) show the extracted IMFs of the synthetic signal (Fig. 1) using the EMD approach and the corresponding FFTs, respectively. The first IMF contains the fastest oscillation content of the data and the last IMF is in fact not an oscillatory signal at all. FFTs show the frequency content moves from the high frequency range to the low-frequency range as the number of IMF increases. This means the lower-level IMFs represent the high-frequency content and the higher-level IMFs represent the low-frequency content. Fig. 5(c) presents the Hilbert spectra of extracted IMFs. In the EMD method, each IMF is ideally expected to be a separate/unique exponentially decaying sinusoid with but this separation may not be realized in real-life signals such as those obtained in SHM. For example, in Fig. 5(c), modes are mixed. If the modes were not mixed the instantaneous frequency plot of the HHT would indicate an almost constant value for each IMF (it is mostly constant for IMFs 4, 5, and 6 but not the others). Another shortcoming of HHT is that it generates spurious modes. For example, the original uncontaminated synthetic signal used in this research has four intrinsic frequencies but the HHT produces 5 oscillatory IMFs.

4.1. Pseudospectrum estimate using Music

Most parametric signal processing methods estimate the signal as a sum of a certain number of complex exponentials. In the methods that represent a real valued signal by a sum of complex exponentials, such as free vibration response of a $d$-DOF system, $2$×$d$ complex exponential components must be found to cover all corresponding dynamics of the system. For example, for any given sinusoidal signal, two complex conjugate exponentials are found whose summation forms the sinusoidal signal. If there is a complex eigenvalue in a real-valued matrix, its conjugate is also an eigenvalue of the matrix. As such, the number of signal components is twice the number of DOFs.

Knowing the number of the signal components, $D=2d\text{,}$ a priori, the multiple signal classification (MUSIC) approach uses estimates of the eigenvectors of a correlation matrix associated with the data vector to calculate a corresponding pseudospectrum [49, 50]. For a given signal of length $N$ (a column vector with $N$ rows) presented as vector $\mathbf{x}$, it is assumed that the signal is defined as the sum of $M$ complex exponential sequences plus noise:

where ${\alpha }_k$, ${\omega }_k$, ${\phi }_k$ are the amplitude, frequency, and phase angle of the $k$th complex exponential sequence and $\mathbf{u}$ is the zero mean white noise signal of length $N$. If $\mathbf{X}$ is a rectangular Toeplitz matrix of $\mathbf{x}$, the covariance matrix of size $N$×$N$ is defined and decomposed as [51, 52]:

where $E$ is the expectation operator used for considering different observations of the signal, $\mathbf{A}$ is the eigenvector matrix of size $N$×$D$, $\mathbf{P}$ is an $D$×$D$ diagonal matrix of the estimates of amplitudes squared, ${\alpha }_k^2$, $\mathbf{I}$ is the identity matrix of size $N$×$N$, ${\sigma }^2$ is the noise variance, and ${\mathbf{A}}^H$ is the Hermitian or conjugate transpose of $\mathbf{A}$ matrix. Then, $\mathbf{x}$ is estimated as:

where vector $\mathbf{p}$ includes estimates of ${\alpha }_k$ with length $M$ such that $\mathbf{P}=E\left(\mathbf{p}{\mathbf{p}}^H\right)$. The eigenvectors corresponding to $M$ largest eigenvalues of ${\mathbf{R}}_x$ span the signal subspace and the eigenvectors corresponding to the remaining eigenvalues of size $J=N-D$ form the noise subspace.

Fig. 5a) Extracted IMFs using EMD method, b) corresponding FFT, and c) Hilbert transforms

a)

b)

c)

The pseudospectrum values are in fact the inverse of the Euclidian distances of all the sinusoidal sequences in the frequency range to the noise space formed by $J$ eigenvectors. The peaks in the pseudospectrum plot versus frequency is expected to reveal the frequency content of the signal. The equation for pseudospectrum as a function of frequency $f$ value is:

where $\mathbf{G}$ is the $N$×$J$ matrix of noise eigenvectors, and $\mathbf{b}\left(f\right)$ is the vector of terms in a geometric progression of exponentials with complex exponent $if$:

To improve the time resolution in this research the pseudospectrum calculation is done for shorter lengths of $\mathbf{x}$ by windowing the main signal with overlapping windows. This way the dimension $N$ is replaced by the window length in Eqs. (10)-(14). Fig. 6 shows the pseudospectrum estimate of the synthetic signal (Fig. 1) using the MUSIC algorithm. The estimated values of the four frequencies are presented in Table 1. The MUSIC algorithm does not estimate the two closely-spaced frequencies accurately.

Fig. 6Pseudospectrum estimate of signal in Fig. 1 using the MUSIC algorithm with the peak frequency values of 0.29, 0.99, 1.57, and 9.47

4.2. Empirical wavelet transform

Empirical wavelet transform (EWT) is an adaptive wavelet transform proposed recently by [53] intended for processing non-stationary signals. In this method the original signal is decomposed into a number of sub-signals representing the modes of original signal utilizing the FT method. The number of wavelet filters and their bandwidths are chosen so that each one works like a band pass filter containing the frequency of one mode of the signal. After the decomposition, each sub-signal is analyzed using the Hilbert spectral analysis to find the time-frequency behavior of each mode.

Fig. 7 shows the FT of the synthetic signal with three detected boundaries through EWT, shown by dashed vertical lines resulting in four frequency regions. The boundaries are detected at the lowest local minima between $d$ (number of DOF) largest maxima. A wavelet filter is generated for each region and is used to create a sub-signal containing a separate mode. In this case, the first boundary in Fig. 7 is not between the first two frequencies resulting in a spurious first detected mode. Further, the region between the first and second boundaries contains the first two modes of the signal resulting in mixing of the first two modes. Fig. 8 shows the four detected components of the synthetic signal and the corresponding Hilbert spectrum using the EWT. The first component is non-oscillatory and represents none of the signal modes because of the incorrect detection of the first boundary. Component 2 represents both closely spaced modes without separating them. Components 3 and 4 represent modes 3 and 4 of the signal but the latter contains high-frequency noise. Performance of this method depends greatly on the correct detection of the boundaries and the scheme used for this purpose is not effective when the signal contains closely-spaced frequencies. EWT yields separate modes in time and time-frequency domains but no system parameter is extracted directly. The estimated values of the four frequencies, the peaks in FT of each output (Fig. 7) are presented in Table 1.

Fig. 7The FT of the synthetic signal using EWT where three boundaries are detected with dashed vertical lines resulting in four frequency regions

Fig. 8a) Four detected components of the synthetic signal and b) the corresponding HT using the EWT

a)

b)

4.3. Approximate Prony method

The idea of representing sequences by sums of damped exponential goes back to mathematician Gaspard de Prony as early as 1795 but was applied in the signal processing practice in the late 1980s [54-56]. Prony’s method was intended for exact representation of noncontaminated signals without any noise [57]. When applied to noisy signals the method exhibits numerical instability and inaccurate frequency and damping estimation and researchers have tried to modify and improve the performance of the method. Beylkin and Monzón [58] used the Prony’s method to efficiently approximate complex exponentials contaminated by noise. The method, called the Approximate Prony method (APM), overcomes the numerical instability of the original method and results in a more efficient approach. Potts and Tasche [59] extended the work of Beylkin and Monzón for non-damped multivariate exponential data and proposed the sparse approximate Prony method that uses a lower number of sampled data. In this research this approach is used for damped exponentials with noise.

The idea of Prony’s method is based on the fact that the solutions of homogeneous linear ordinary differential equations are sum of complex exponentials. In addition, linear differential or difference equations are transformed to a linear recurrence sequence, similar to autoregressive (AR) models. Therefore, the signal must fit in a recurrence model. The roots of the polynomial formed by the coefficients of the AR model describe the model in explicit form. In Prony’s method the signal is assumed to be sum of complex exponentials, and is decomposed into a preselected number of exponentially decaying sinusoids in the following form:

where:

${\alpha }_n$, ${\mu }_n$, $f_n$, ${\phi }_n$ are the amplitude, damping exponent, frequency, and phase angle of the $n$th component respectively, $S_n$ is called the residue, weight, phasor (in electronics) or complex coefficient, and $z_n$ is called the pole corresponding to the $n$th component of the signal since they are the roots of a polynomial to be described in the next paragraph.

Using the number of signal components, $D$, and the signal length $N$, an integer $L$ is chosen in the range $D<L<N/2$ depending on the types of the signal and its noise content and a data Hankel matrix (a flip of Toeplitz matrix), ${\mathbf{F}}_1$, and a data vector, $\mathbf{v}$, are formed similar to AR model of order $L$ as follows:

Next, the following equation is solved to find the constant complex coefficients vector $\mathit{c}={\left[c_L,\dots ,c_1\right]}^T$:

Therefore, assuming $N>2L$:

where ${\mathbf{F}}_1^{+}$ is the Moore-Penrose pseudoinverse of ${\mathbf{F}}_1$ used to solve the linear least squares problem defined by Eq. (20). Coefficients of the recurrence sequence, $c_i$, are used to form an $L$th order complex polynomial as follows:

Roots of Eq. (22) are the estimates of poles defined in Eq. (17). In noiseless situations, only $D$ non-zero coefficients are found by Eq. (21) that completely describe the signal content. However, with the existence of noise, $D$ roots are chosen that best match the signal components among possibly more than $L$ non-zero roots. Next, using the $D$ selected roots, the $N-L-1$ by $D$ Vandermonde data matrix or the matrix of found exponential columns is formed and used to find the complex coefficients (residues) vector $\mathbf{S}$ to be used in Eq. (15):

The main challenge in this method is the selection of the roots that are related to signal components and not spurious modes.

Table 1 presents the estimated values of the four frequencies and their corresponding damping exponents. This method tends to distinguish the two closely-spaced frequencies but not accurately at 0.2 and 0.55 Hz compared with the actual values of 0.3 and 0.5 Hz. The inaccuracy is reflected in the estimated signal shown in Fig. 9 which does not replicate the original signal accurately. Further, there exists a spurious frequency at 9.81 Hz that contributes to the high-frequency oscillation during the first second of the estimated signal, with a high damping exponent that contributes to the fast dissipation of this mode as seen in Fig. 9.

Fig. 9Comparison of the signal estimated by the AP method with the original synthetic signal

4.4. Matrix Pencil method

Similar to the AP method, the Matrix Pencil (MP) method is also used to decompose the signal into a preselected number of exponentially decaying sinusoids. For any given signal, an integer $L$ called the pencil parameter is chosen in the range $D<L<N/2$ and two data matrices ${\mathbf{F}}_1$ and ${\mathbf{F}}_2$ of size $(N-L)\times L$ are formed similar to data matrices in the Prony’s method but with a delay of one time increment as follows [60, 61]:

Parameter $L$ is chosen to minimize the effect of noise in the data as explained in the AP method. Assuming linear approximation, matrices ${\mathbf{F}}_1$ and ${\mathbf{F}}_2$ are written in the following form:

where:

where ${\mathbf{Z}}_L$ and ${\mathit{Z}}_R$ are $(N-L)\times D$ and $D\times L$ matrices, respectively, and $S_i$ and $z_i$ are defined by Eqs. (16) and (17). The Pencil of matrices ${\mathbf{F}}_1$ and ${\mathbf{F}}_2$ is:

Eq. (32) shows that if $\lambda$ is equal to any of poles ($z_i$), it is also an eigenvalue of the pencil of matrices ${\mathbf{F}}_1$, ${\mathbf{F}}_2$. As such, finding eigenvalues is equivalent to finding the poles of the signal as:

where superscript “$+$” denotes the pseudoinverse of the rectangular matrix defined as:

Table 1 presents the estimated values of the four frequencies and their corresponding damping exponents. This method distinguishes the two closely-spaced frequencies accurately at 0.32 Hz and 0.55 Hz compared with the actual values of 0.3 Hz and 0.5 Hz and generates no spurious solutions. Fig. 10 shows comparison of the signal estimated by the MP method with the original synthetic signal. The MP method results in a better match in time domain compared with the AP method confirming the earlier assertion by Sarkar and Pereira 1995 that MP works better in the presence of noise without actually comparing the two methods.

Fig. 10Comparison of the signal estimated by the MP method with the original synthetic signal

4.5. Esprit method

Estimation of Signal Parameters by Rotational Invariance Technique (ESPRIT) [62] explores the rotational invariance property in the signal subspace. The guiding principle of ESPRIT is closely related to that of Matrix Pencil [63, 64]. ESPRIT estimates the frequencies ${\omega }_n$$(i=1,\dots ,D)$ as $\arg \left(v_i\right)$ (for argument or phase angle) where $v_i$s are the eigenvalues of the matrix $\mathrm{\Phi }$ as follows:

in which $S_1$ and $S_2$ are data matrices similar to $F_1$ and $F_2$ defined by Eqs. (29) and (30).

The frequencies estimated by this approach for the synthetic signal are presented in Table 1. ESPRIT is also unable to detect the closely spaced frequencies.