Walsh Transform and Empirical Mode Decomposition Applied to Reconstruction of Velocity and Displacement from Seismic Acceleration Measurement

This paper focuses on reconstruction of dynamic velocity and displacement from seismic acceleration signal. For conventional time-domain approaches or frequency-domain approaches, due to initial values and non-negligible noise in the acceleration signal, drift and deviation in velocity and displacement are inevitable. To deal with this deficiency, this paper develops a Walsh transform and Empirical Mode Decomposition (EMD)-based integral algorithm, or WATEBI in short. In the WATEBI algorithm, the Walsh transform is employed to realize vibration signal reconstruction. Next, the EMD method is used to eliminate the residual in the reconstructed signal. Finally, the trend term in velocity and displacement is removed by linear least-squares fit. This algorithm can be straightforwardly implemented by an ordinary computer. Reconstructed displacements and velocities from vibration of a simulated single-degree-of-freedom system and two-site measured ground motions in earthquakes validated the robustness and adaptiveness of this algorithm. It can be also applied to many other areas, like mechanical engineering and ocean engineering.


Introduction
With the rapid development of modern science and technology, the research on mechanical and structural vibration has promoted great advances in numerous engineering fields and areas, such as aerospace engineering, bridge engineering, earthquake engineering and offshore engineering [1,2]. In mechanical engineering, vibration tests have attracted considerable attention, as they are of significant value for the health monitoring of machinery [3], including system operational reliability and mechanical breakdown detection. In earthquake engineering, ground motion signals are the key to evaluating seismic magnitude and intensity [4]. Records of ground motions in a strong seismic event are an essential input to time-domain dynamic response calculation for structures located in earthquake-prone zones. Thus, it is important to obtain accurate vibration responses, either directly through in situ real-time monitoring, or indirectly through vibration tests and numerical reconstruction through data post-processing. Usually the responses of interest are acceleration, velocity, displacement, and stress signals of a structure. They can be measured by a wide variety of sensors. In most situations it is acceleration that is measured. Please note that the measured acceleration is an absolute quantity, while the measured displacement is relative to a reference object. In reality, it is not straightforward to install a fixed displacement sensor due to the geometrical complexity of a structure and its restricted ambient space [5]. By contrast, as a result of their small size, acceleration sensors need neither a large space nor a stationary reference [1]. Among concurrent motion sensors, acceleration sensors (for instance accelerometers) are technically the most versatile. It is very convenient in both laboratory tests and field observations to carry out acceleration logging.
To acquire data of unknown velocity and displacement, a traditional and popular approach is the integral of measured acceleration. There are two main integral schemes. One is by integrating the circuit and the other uses a numerical integral algorithm [6]. Circuit integration is rarely used as it is sensitive to the low frequency of acceleration, while the numerical integral scheme has been widely adopted for years. The numerical scheme can be divided into two categories, time-domain integral [5] and frequency-domain conversion [7]. The former includes known methods like trapezoid formula, Simpson's rule, Newton-Cotes formula and first to Nth-order integral operator based on the Runge-Kutta algorithm. However, these conventional methods have several disadvantages [1]. The result usually exhibits drift and errors in the trend term [8] because of the contained direct current component, and high-frequency or low-frequency noise. The frequency-domain conversion method by contrast is numerically much faster, by transforming the time-domain acceleration into a curve synthesizing complex coefficients through finite Fourier series. Various Fast Fourier transform (FFT) algorithms have been widely applied to build the corresponding relationship between frequencies and amplitudes [9]. To solve the problem of spectral leakage in FFT, the antileakage least-squares spectral analysis (ALLSSA) [10] and multichannel ALLSSA [11] have been proposed to regularize the irregularly sampled series and to attenuate the random noise. They were proven more robust than the antileakage Fourier transform (ALFT) [12] and the arbitrarily sampled Fourier transform (ASFT) [13]. It follows that by the principle of trigonometric function integrals, the velocity and displacement signals in frequency domain could be obtained by simply dividing the Fourier coefficients of acceleration with transfer functions (iω) and (iω) 2 . Then, the time series of velocity and displacement are recovered by applying the Inverse Fourier Transform (IFT) [7].
Many attempts have been proposed to improve the rationality and accuracy of the aforementioned methods. These attempts are categorized into three classes of integral approaches, the time-domain polynomial fitting to remove the trend error [14], the band-pass filtering method [15], and FFT implementation with low-frequency cutoff and attenuation attempts [16][17][18], have demonstrated their good applicability in engineering problems. Nevertheless, these approaches still have flaws with respect to error control and fidelity. Recent years have also seen the application of Hilbert-Huang transform and wavelet transform [19][20][21][22], such as the least-squares wavelet and crosswavelet analyses [23]. However, these transforms have not yet been well established for signal conversion.
This paper proposes a Walsh Transform and Empirical Model Decomposition (EMD)-based Integral algorithm (WATEBI) to obtain the desired velocity and displacement from acceleration. To demonstrate the level of robustness and adaptiveness of this algorithm, three case studies are presented. The first study is a numerically simulated vibration of a single-degree-of-freedom (SDOF) system driven by white noise. The other two cases are seismic recordings including ground acceleration, velocity and displacement, respectively, measured in two typical earthquakes.

Walsh Series and Transform
The Walsh series is a family of non-sinusoidal orthogonal rectangular functions, originally defined by Walsh in 1923 [24,25]. Figure 1 shows the first eight Walsh series. Each function in the Walsh series, of different orders, jumps from 1 to −1 within the interval [0, 1] [26]. They are orthogonal, symmetrical and linear, thus extremely suitable for signal reconstruction [27], but only a few studies in literature have addressed their applicability to the recovery of random series. In digital processing, the discrete Walsh transform (DWT) and the inverse DWT (IDWT) of a signal x(t) in terms of the Walsh series wal(k, n) are defined as [28]:

Walsh Transform versus Fourier Transform
Unlike the Fourier transform, the Walsh transform is operated by real number additions and subtractions. Walsh functions are formed in rectangular waveforms of ±1 that are piecewise constant. This characteristic makes the calculations much less complicated, as opposed to complex number handling in FFT [29,30]. Moreover, all Walsh functions are arranged by a certain order and in constant formations ( Figure 1). Their shapes are invariable regardless of the selected order number k for a given random signal. As a result, the stability of the integral result is guaranteed. On the contrary, the Fourier transform aims to convert a signal into a trigonometric series whose frequencies, amplitudes and corresponding resolution depend on the sample number N used in FFT. In addition, in accordance with the Nyquist-Shannon sampling theorem [31], the Nyquist folding frequency fc in FFT is determined by the time interval Δt. Aliasing of components higher than fc into the spectral zone of f < fc is inevitable. Hence, selections of parameters would significantly influence the integral results [17].

Empirical Mode Decomposition (EMD)
The Empirical Mode Decomposition (EMD) method is adopted in the WATEBI algorithm for removing the nonlinear noise in the reconstructed velocity and displacement. The EMD method was proposed by Huang et al. [32]. The vibration signal can be decomposed into several intrinsic mode functions (IMF) and a residual term shown as Equation (3). where

IMF[x(t)] are the intrinsic mode functions and RS[xr(t)] is the residual term.
It is suggested that the remnant RS[xr(t)] of a signal serving as random inutile residue be discarded after the Hilbert-Huang transform [32][33][34]. Thus, the last term of EMD is removed from the velocity and displacement signal in the WATEBI method.

Signal Reconstruction
A digital time history of acceleration can be first transformed into N-order Walsh series as follows.
where a(n) is the varying acceleration at time instant of (n-1)Δt, Δt is the equidistant time interval. The velocity is then obtained by integration of the right-hand-side (RHS) once, despite that its initial condition might be unknown. where is the discrete version of the integral of Walsh series [27]. Figure 2 illustrates the first four integrals of Walsh series. Obviously, they are triangles and linear.
A(k) and V(k) in Equations (4) and (6) are respective Walsh coefficients of acceleration and velocity, indicating the amplitudes of Walsh series in Figure 1. These two equations both involve integrals of wal(k, n) that are identical and therefore need to be computed only once. A 2-dimensional N-by-N matrix can be used to store them beforehand. Apparently, the above algorithm is numerically simple and fast as it only requires superposition of the integrals of series wal(k, n) with respective coefficients for k = 0, 1, …, N − 1 that can be easily fulfilled by matrix manipulations. In this study, the MATLAB-based Hadamard-Walsh transform [35] is used to obtain the Walsh coefficients. The residual noise in the discrete velocity and displacement is then removed by the EMD method respectively. Figure 3 illustrates the procedures of implementing WATEBI algorithm.

Numerical Effectiveness and Elimination of Trend Error
It should be noted that the amplitude, namely coefficient A(k) or V(k), of a high-order Walsh series is very small, such that the corresponding integral contributes little to the overall displacement. Additionally, analogous to Fourier series, high-frequency components in the original signal are usually regarded as noise that ought to be filtered off. Therefore, in the WATEBI algorithm, only the low-order series are retained for use. By doing so, the computational time is reduced by approximately half as well. In general, the first half of Walsh series suffices the accuracy requirement while the second half are truncated, i.e.
Considering that this WATEBI algorithm only uses first-order integrals of Walsh series in the linear formation, the accumulation of integral errors of the derived velocity is insignificant.
For displacement, although the linear superposition of Walsh series is performed for one more time, the integral kernels in Figure 2 are unchanged. Hence, the accumulated error is one order lower than the square term from conventional time-domain integration methods. Nevertheless, the trend error still exists and should be removed in order to eliminate its subsequent influence on the displacement. In this study, the least-squares method is employed to fit the linear term in trend error of velocity and displacement.

Case Studies and Analyses
To assess the performance of the proposed WATEBI method, three case studies are presented in this paper. These are vibration of a SDOF system and the ground motions of El Centro earthquake and Kaikoura earthquake. In each case study, the WATEBI method is compared with time-domain trapezoidal integral by Labview [36], MATLAB Ode45 (Ordinary Differential Equation) function with tolerance of 10 −6 [37], low-frequency cutoff [18], and low-frequency attenuation with the target accuracy αT = 0.97 [16,17]. For a fair comparison among all methods, both the linear trend in velocity and the nonlinear trend in displacement have been removed, respectively by linear fitting and nonlinear fitting. The only difference is that the WATEBI method employs EMD while the benchmark methods use the common least-squares method.
Please note that the attenuation algorithm in reference [16] is only given for displacement. The accuracy of the frequency domain methods (low-frequency cutoff [18] and low-frequency attenuation [16,17]) heavily depends upon parameter selection. In this paper, the parameters in the selected frequency domain methods for all three case studies are tuned to obtain the ideal results. In Case 1, the cut-off frequency in the cutoff method and the target frequency in the attenuation method are 0.1 Hz, coincidently. In Case 2, they are 0.07 Hz and 0.1 Hz respectively. In Case 3, they coincide again at 0.15 Hz.
In each case study, the following three indicators are used to evaluate the deviation of the reconstructed signal away from the reference velocity/displacement. (13) Average peak error Er1, which denotes the average deviation at both positive and negative extremes of the time series, is defined as where x(t) is the reference signal and X(t) is the reconstructed signal.
(2) Maximum peak error Er2, denoting the larger deviation at positive and negative extremes of the time series, which is often the most concerned and useful information for comparison.
(3) Root-mean-square error Er3 represents the overall deviation of the whole time series at every time instants.

Case 1: Vibration of a SDOF System
A SDOF vibration system driven by white noise is considered first.

⋅ ( ) + ⋅ ( ) + ⋅ ( ) = − ( )
where x, v and a are respectively displacement, velocity and acceleration of system vibration. The assumed frequency of excitation F ranges from 0.1 Hz to 5.1 Hz and its power spectrum is constant, 1 N 2 s. Mean of F is zero. The SDOF system has a vibration mass M of 1 kg, natural frequency fn of 0.45 Hz and damping ratio of 3%. The initial velocity and displacement are assumed as 0 cm/s and −3.2 cm. The dynamic responses including acceleration, velocity and displacement are solved by the Newmark-β method of β = ¼ [38]. This method is stepwise and has been widely applied in civil engineering for years due to its accuracy, as long as a suitable interval Δt is chosen in computation. This study uses a small time interval, Δt = 0.0586 s. Figure 4 shows the time-domain numerical solution of excited acceleration for duration of 239.94 s.  Figure 5a shows the velocities reconstructed by the five methods against the reference velocity solved by the Newmark-β method. The comparison of displacements is given in Figure 6a. In these two figures, it can be clearly seen that after removal of trend by fittings, none of the reconstructed signals contains drift. Moreover, at the beginning and ending parts of the full duration, these two frequency domain methods fail to precisely reconstruct both velocity and displacement. By contrast, as illustrated in the enlarged details in Figures 5b and Figure 6b, the velocity and displacement reconstructed by the WATEBI method match the reference time series with smallest deviations. The comparisons of average peak, maximum peak error and root-mean-square error for the underlying five methods are given in Tables 1 and 2, respectively for velocity and displacement. It can be shown from the three error indicators that for velocity the low-frequency cutoff method performs worst, while for displacement the low-frequency attenuation is of lowest accuracy. For both velocity and displacement, the best performance among all five methods is WATEBI. In addition, the good performances of Matlab Ode45 and trapezoidal integral are noticeable, because they are analogous to the Newmark-β algorithm for the adopted step-by-step integration scheme in the time domain and actually in this SDOF system no noise is contained in the acceleration signal. Nonetheless, in the following case studies, the limitation and drawback of time-domain methods will be demonstrated when random noise is mixed in the acceleration.

Case 2: El Centro Earthquake Ground Motion
The El Centro earthquake ground motion, widely used in civil and earthquake engineering, was obtained from the Imperial Valley earthquake in 1940, recorded at the site of El Centro. The acceleration record is a 2688-point signal with Δt = 0.02s. It can be seen from Figure 7a that the maximum ground motion exceeds 300 Gal (cm/s 2 ). Figure 7b shows the amplitude spectrum of acceleration. Obviously, it is wide-banded, and contains numerous frequency components as many other seismic records do and the frequencies of significant components are lower than 10 Hz. Noise may exist in the record due to low-frequency tilt of device, high frequency burr or even direct component, all possibly leading to deviation in the reconstructed velocity and displacement.
The site-recorded velocity and displacement are used as reference for assessing the accuracy of WATEBI algorithm. Before running WATEBI, the original acceleration history is interpolated with the cubic spline function using Δt/6, in order to achieve a finer and smoother curve. Although the original ground motion of acceleration is stochastic and transient, it can be seen from Figures 8 and 9 that the reconstructed signals (especially velocity) by WATEBI approach more closely to the site records than the other four methods do. For displacement, the deviation from the real data starts at some 35 s, but the two most useful peak-to-trough cycles (from 12 s to 28 s) for checking structural dynamic responses are well matched. Tables 3 and 4 compares the performance of five methods for numerically derived velocity and displacement. It is obvious that the trapezoidal integral is still efficient in velocity reconstruction, with regard to error control on peaks, but for the overall performance of velocity and displacement, the proposed WATEBI method outruns all the rest four methods. It should be also recognized that the three error evaluating indicators (Equations (10)-(12)) of displacement results are larger than those of the velocity results. After all, displacement reconstruction is based on velocity. Both noise and the tiny trend term in velocity can amplify the deviation in displacement.

Case 3: Kaikoura Earthquake in The South Island of New Zealand
The Kaikoura earthquake that occurred in the South Island of New Zealand on 13 November 2016 was a very strong seismic event with magnitude of 7.8 (Mw) and a long duration of approximately two minutes. The seismic data used in this paper is a horizontal component that has been corrected after occurrence. The data was recorded with a constant time interval of 0.005 s. Figure  10 shows the time history of acceleration and its amplitude spectrum. Compared with the El Centro record in Figure 7, the energy of Kaikoura record is more concentrated, not only for time history but also for spectrum. Owing to the fine time interval, the total number of logging is 14,327, such that interpolation of data is unnecessary. The reconstructions by WATEBI only require seconds of computational time, using an ordinary computer. Figures 11 and 12 illustrate the waveforms of original and derived velocity and displacement. The good performance of WATEBI algorithm is again demonstrated. The low-frequency cutoff method shows the worst degree of reconstruction among all methods, especially in displacement. Considering that more and more advanced technologies have been applied in earthquake monitoring, rich and accurate seismic recordings are available. This in turn helps to improve the accuracy in deriving velocity and displacement in terms of acceleration. Therefore, very good agreement between the derived and recorded data can be expected.  Tables 5 and 6 give the comparison of five methods in deriving the velocity and displacement of Kaikoura earthquake. Again, the three error indicators show that the proposed WATEBI algorithm is most robust. For displacement, both the average peak error and the maximum peak error by WATEBI are even below 1%. The low-frequency attenuation method also provides good accuracy, because the target frequency has been tuned according to the target accuracy T = 0.97. By contrast, the low-frequency cutoff is the worst among all five methods, not only for velocity but also displacement.

Concluding Remarks
Accelerometers are usually employed to measure vibration signals in engineering practice. To numerically derive random velocity and displacement from recorded acceleration, a Walsh Transform and Empirical Mode Decomposition (EMD)-based integral algorithm, WATEBI in short, is proposed in this study. It is actually the linear superposition of the integral kernels of Walsh series and can be handled by an ordinary computer with ease. The residual noise in the reconstructed vibration signals is eliminated by EMD. Unlike the popularly used low-frequency cutoff and lowfrequency attenuation methods, this algorithm does not need to specify the critical parameters like cutoff frequency, target frequency and target accuracy. As shown in the flowchart (Figure 3), the trend term is removed by linear least-squares fit.
To test the performance of the WATEBI algorithm, three distinct case studies are used. They are respectively the vibration of a SDOF system with given parameters, ground motions of El Centro earthquake and Kaikoura earthquake. In each case study, all reference velocity and displacement are available. Then for the comparison among five methods, trapezoidal integral, Matlab Ode45, lowfrequency cutoff, low-frequency attenuation and WATEBI, the reconstructed signals of velocity and displacement are investigated. The average peak error, the maximum peak error and the root-meansquare error, are three criterions used to differentiate the level of quality of these methods. The graphical comparisons of both full duration and enlarged details in Figures 5, 6, 8, 9, 11 and 12 reveal that not only the time history of reconstructed velocity by WATEBI perfectly overlaps the reference signal, but also the reconstructed displacement by WATEBI tallies with the reference data with very little deviations, especially at peaks and troughs. The errors in Tables 1-6 also demonstrate that for both velocity and displacement, the overall performance of WATEBI is the best among all compared methods. Particularly for displacement, WATEBI behaves very stably, providing lowest errors for all three cases. The Ode45 algorithm works very well in Case 1 and Case 3 for velocity reconstruction, while the trapezoidal integral method is suitable for all cases. The low-frequency cutoff method gives accurate velocity only in Case 1, while the low-frequency attenuation method renders accurate displacement in Case 3. Nevertheless, none of these four methods surpass the WATEBI algorithm in nearly all cases. Additionally, the WATEBI algorithm is more adaptive, since, as depicted in the flowchart (Figure 3), the signal reconstruction does not require subjective parameters such as cutoff and target frequencies to be specified. Signals adopted in these cases are equally spaced. However, for the irregularly sampled data series, the ALLSSA method might be more robust for the preprocessing than higher-order spline methods.
The merit of Walsh transform with constant kernels has inspired this research work. For the sake of limited length of this paper, many other case studies to validate the robustness of the proposed algorithm are not included. Different from the long and popularly used Fourier transform and Hilbert transform in signal processing, the Walsh transform still looks like a stranger to many researchers. Hopefully more virtues of the Walsh transform could be rediscovered to further invoke its wider application in engineering problems.