Automated Denoising Technique for Random Input Signals using Empirical Mode Decomposition (EMD)-Stabilization Diagram

The present paper deals with the novel approach of filtering technique using hybrid of empirical mode decomposition technique with stabilization diagram, that autonomously implemented within Matlab. Noise or unwanted signal is always present in the data and a bad signal-to-noise can cause a severe error in modal parameter extraction. With the recent developments of automated procedures without user interaction for the operational modal analysis (OMA), the corrupted input signals turn out to be a big issue in obtaining reliable results of automated modal parameter identification. The appearance of noise or unwanted modes due to environmental effects could affect the actual structural modes selection. There is a significant issue regarding “noise” (or spurious) modes and eliminating them from the raw signal remains to be solved and requires a lot of interaction with an expert user. In the parametric modal analysis, oversizing of a modal model is usually performed to minimize the bias on modal estimates by getting all physical modes in the frequency range of interest and help to obtain a good model fit to the data. However, this will introduce noise modes. Thus, authors take advantage of tools, such as the stabilization diagram, to illustrate, and decide, if a mode is physical or not. This selection is not a trivial task, but it may be difficult and time consuming depending on the quality of data, the performance of the estimator and the experience of the user. Since the extensive interaction between tools and user is inappropriate for monitoring purposes, image clustering tool is introduced to separate the physical poles from the others with short response time and low computational efforts compared to the available clustering algorithm. Meanwhile, Empirical mode decomposition (EMD) is then introduced to break down a signal into various components without leaving the time domain and purposely used for filtering. These are a great combination as well as an effective procedure in producing a good input signal that free from unwanted modes which can cause disruptive decision making for the actual modes selection.


Introduction
Noise or unwanted signals are always present in the data and a bad signal-to-noise can cause a severe error in modal parameter extraction. With the recent developments of automated procedures without user interaction for the operational modal analysis (OMA), the corrupted input signals turn out to be a big issue in obtaining reliable results of automated modal parameter identification. The appearance of noise or unwanted modes due to environmental effects could affect the selection of the actual structural modes.
Several examples of filtering technique performed in OMA have been reported in the technical literature and presented herein. Yang et al. proposed using a band-pass filter based on Hilbert-Huang transform (HHT) in modal identification method via measured free vibration time histories [1,2]. Yu and Ren proposed an EMD based SSI modal identification method, which reduced the postprocessing work in SSI algorithm (picking the stable poles) and using a band-pass filter to treat the closely spaced modes [3]. Among the methods based on conventional signal processing, the so-called time domain filtering method was a tracking procedure based on the application of a band-pass filter to the system response in order to separate the individual modes in the spectrum [4].
Abovementioned literature indicates that the researchers often used classical filtering technique based on Fourier transforms, particularly band-pass filter rather than low pass and high pass filter in order to remove unwanted signal and also for clear modal parameter extraction. It is possible to remove unwanted from its raw signal, but more attention should be paid to ensure all the necessary information from this signal are not affected. Therefore, it is still worthy to investigate this problem.
In the parametric modal analysis, oversizing of a modal model is usually performed to minimize the bias on modal estimates by getting all physical modes in the frequency range of interest and help to obtain a good model fit to the data [5]. However, this will introduce noise modes. There is a significant issue regarding "noise" (or spurious) modes and eliminating them from the raw signal remains to be solved and requires a lot of interaction with an expert user.
Thus, this paper will focus on automated OMA that can ensure an effective identification of physical modes and removing unwanted signal with short response time and low computational efforts.

Automated filtering technique
The automatic filtering technique in OMA based on parametric methods consists of the following steps: (1) Measure the responses of the structure and estimate the modal parameters using a high model order, n modes.
(2) Construct Stabilization diagram by estimating poles with an increasing model order and illustrate, and decide, either mode is physical or computational modes (3) Classify the n modes in physical and computational modes using a clustering algorithm. (4) The unstable mode is filtered out involving empirical mode decomposition (EMD) method and left only the useful signal.
The acceleration time series of the first data set from the Heritage Court Tower building has been processed by the automatic versions of the SSI-COV. The frequency range of 0 to 5 Hz already contains a considerable number of resonant frequencies.

Stochastic subspace technique (SSI)
The Stochastic subspace identification (SSI) has been well-known over the past decade, mainly due to userfriendly executions [6]. There are two main types of SSI, that is, the correlation driven, COV-SSI and the data-driven, DATA-SSI. The DATA-SSI is an approach where the estimation of the correlation functions is integrated into the identification algorithm. The COV-SSI, however, is very alike to the ERA [7] and identifies a stochastic state space model from output-only data [5]. This paper only concerns on correlation driven SSI and details derivation about this method are described below.
The first step of this method is to calculate the output correlations (Equation 1). denotes the unbiased estimate of the correlation matrix at time lag i based on a finite number of data: (1) Where is the data matrix Y with the last block rows i removed and is the transpose data matrix with the first block rows i removed. Therefore each matrix gets dimensions l*l. The estimated correlations at different time lags (Equation 2) are then gathered into a matrix called the block Toeplitz matrix: (2) The Toeplitz matrix contains i*i matrices and is therefore of dimensions li*li. For the identification of the modal parameters of a system of order n, the Toeplitz matrix needs to be n*n. Therefore, the following need to be true for the number of block rows i: (3) Where n is the system order. Normally, the system order is unknown and the following (Equation 3) should be fulfilled. For complicated structures the number of block rows, i should be higher than this minimum criterion for better results. The magnitude x of this relies on the problem and should be chosen as an input by the user. In this case, the magnitude of block rows, x was set equal to 2, meanwhile, for maximum system order, n was defined as 50 modes.
This value needs need to be selected wisely. If the maximum order is chosen to be smaller than the correct system order you will not get any correct results. If you, however, assume this value too high you will get too many nonphysical modes and it will be harder to derive which modes that are the actual correct physical modes, in addition, the computational time greatly increases.
Then the singular value decomposition (SVD) of the block Toeplitz is calculated. It gives the unitary  [8]. (4) To extract the dynamic response, you need to obtain the state matrix [A]. This is done for each order from 1 to . The first thing needed is to find the observability matrix and the reversed controllability matrix which is found by the factorization of .
Where the observability matrix with dimensions li*n is given by: And the reversed controllability matrix with dimensions n*li is given by: The SVD of is already computed in Equation 4 and we can use the result to find and . By splitting the SVD into two parts and using the identity matrix [I]: The mode shapes of the system are obtained from and : When the modal parameters are found for the current order, the process of finding the state matrix and the modal parameters are repeated for each order up until . Then all the parameters are plotted in a stabilization diagram, which will be discussed in the next section.

Stabilization Diagram
When the modal parameters are attained, the stabilization diagram can be constructed. This tool is one of the most common methods to distinguish physical poles from others, which is performed by estimating poles with an increasing model order. Since the model of the system is often oversized, thus, the plot will contain noise modes and mathematical modes. The noise modes are caused by physical reasons, while the mathematical modes are generated to ensure the mathematical description of the measured data. Theoretically, the physical poles should stabilize and can be found by an alignment of stable poles, whereas the computational or mathematical poles are scattered, showing the criterion of the unstable poles. This is based on the comparison of the poles associated with a given model order with those attained from a one-order lower model [5].
The natural frequencies and damping ratio of poles from two orders are compared [10]: Only the poles that fulfill a stabilization criterion defined by the user (x and y) are labeled as stable respected the following limits for variations between models of consecutive orders: natural frequency variation < 1% [11]; modal damping ratio variation < 5%. These limits proved to permit the clear distinction of vertical alignments of stable modes that represent the modes of vibration.

Image Clustering
The poles of the stabilization diagrams that represent the same physical modes were easily grouped by using image clustering. The process of image clustering will be explained in this section.

Input image
The stabilization diagram was cut and displayed separately into every frequency according to 0.01 interval, 5Hz / 0.01 = 500 images. So, every image represents the frequency of 0.01 Hz.

Image feature extraction
In this study, six standardized image features in Matlab software were used to extract the image features of each stabilization diagrams. These features specifically represent the characteristics of each parameter (natural frequencies, damping ratios) for different conditions either stable or unstable. Table 1 summarises all the image features that were extracted in this study and their characteristic values. This algorithm was configured to group all the poles of the stabilization diagram, stable and not stable, by using image features extraction, as shown in Figure 3.

Empirical mode decomposition (EMD) method
Empirical mode decomposition (EMD) is an effective approach to decompose a nonstationary signal into various components, known as a series of intrinsic mode functions (IMFs) and a residual, without leaving the time domain [12]. The basic steps of EMD are as follows.
a) Sifting process: Find all the extremes of ( ). Define the upper and lower envelope of ( ), which can be obtained through connecting all the local maxima and local minima with cubic a spline curve, respectively. Calculate the mean value of the two envelopes. The difference between the ( ) and the mean value is designated as ℎ1( ). b) Checking process: Check whether ℎ1( ) satisfies the two conditions of IMF. If ℎ1( ) is an IMF, then take 1( ) = ℎ1( ) as the first IMF of ( ). Otherwise, let ( ) = ℎ1( ) and repeat the sifting process until the first IMF is found out. c) Looping process: Define 1( ) = ( ) − 1( ) as the new signal and repeat (a), (b) until ( ) is smaller than a predetermined threshold. Then, the original signal ( ) can be rewritten as The acceleration time series signal of the first data set from the Heritage Court Tower building has been breaking down into 13 IMF as shown in Figure 4 below. This figure shows the all the IMFs of the acceleration response signal. The response signals were decomposed into time-domain and of the same length as the original signal allows for varying frequency in time to be preserved, particularly in decreasing order.

Results and discussion
The selection of the weak modes of the system that autonomously implemented involved Matlab command -find and the threshold in order to discriminate the actual modes and left only the weak ones. The approach is based on excluding the dominant modes with a higher value than the threshold. The threshold is determined by the half of maximum peak in image clustering plot. Meanwhile, for the selection of useful IMF and unwanted IMFS are based on the frequency of unwanted mode. The peak of a certain frequency in IMFs can be determined by using peak picking (PP) method. Thus, the useful IMF will reconstruct back to become an input signal. The result of the automated filtering technique using image clustering and EMD is characterized in Figure 5. From the image clustering plot above we can see a peak at frequency 0.21 Hz had been significantly reached a low point after went through filtering process when compared to image clustering plot in Figure 3, where a shoot up peak clearly seemed at a similar frequency. Thus, the results of this study indicate that, the implementation of image clustering for identification physical modes and computational modes in stabilization diagram and EMD as a tool for filtering purpose are a great combination as well as an effective procedure in producing a good input signal that free from unwanted modes which can cause disruptive decision making for the actual modes selection. The use of image clustering rather than normal clustering algorithm was capable to reduce processing time and low computational efforts needed. In the case of application of advanced clustering algorithm, the number of parameters had to be set at startup which consumes time for calibrations in order to allow reliable identification of structural modes. Meanwhile, EMD is a suitable method for filtering based on the frequency of peak because it decomposes a signal into various components without leaving the time domain and according to decrease order of frequency. Thus, the unwanted frequency peak can easily be removed out from the original input signal. However, this method has drawbacks in separating closely space modes.
This paper has given an account of and the reasons for the widespread use of stabilization diagram, image clustering and EMD for automated filtering technique that involving random input signal. These findings provide the following insights for future research in term of eliminating noise or unwanted modes from closely spaced modes. Even though, image clustering was capable to identify closely spaced modes, as appeared in Figure 1 for the all closely spaced modes of the Heritage court building from 1.2 to 1.5 Hz. But, more attention should be paid to select the appropriate tool or improving existing EMD in eliminating this closely spaced noise or unwanted modes.