# An Underdetermined Blind Source Separation Method with Application to Modal Identification

Authors: Gang Yu
Journal: Shock and Vibration
 (2019)

---

## Abstract

In structural dynamic analysis, the blind source separation (BSS) technique has been accepted as one of the most effective ways for modal identification, in which how to extract the modal parameters using very limited sensors is a highly challenging task in this field. In this paper, we first review the drawbacks of the conventional BSS methods and then propose a novel underdetermined BSS method for addressing the modal identification with limited sensors. The proposed method is established on the clustering features of time-frequency (TF) transform of modal response signals. This study finds that the TF energy belonging to different monotone modals can cluster into distinct straight lines. Meanwhile, we provide the detailed theorem to explain the clustering features. Moreover, the TF coefficients of each modal are employed to reconstruct all monotone signals, which can benefit to individually identify the modal parameters. In experimental validations, two experimental validations demonstrate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

One of the main issues in structural dynamic analysis is to identify the modal parameters, e.g., frequency, damping, and modal shape. Starting from [1, 2], the blind source separation (BSS) technique has become more and more popular in modal identification, due to its straightforward, computationally fairly efficient, nonparametric, and requiring no prior information of the dynamic system. From early literatures, it can be known that the applications of conventional BSS methods, e.g., independent component analysis, second-order blind identification, and their improved versions, mainly focus on fundamental research and theorem analysis [3–7]. As research interest increased, some limitations of conventional BSS techniques are gradually recognized. An essential assumption that guarantees successful application of conventional BSS methods is that the number of sensors should be not less than that of sources, namely, the issue of the determined BSS. However, considering various limitations in practice, e.g., costs, and difficulties in access, the installation of multiple sensors may not be feasible. Furthermore, without the prior knowledge of active modals in most cases, the requirement of adequate sensors is too strict to meet [8, 9]. Therefore, there is a strong requirement to develop the underdetermined BSS techniques for structural dynamic analysis, which can address more sources with limited sensors.Recently, an underdetermined BSS method called sparse component analysis (SCA) draws many attentions in structural dynamics [10, 11]. By transforming the original signal from time domain into sparse domain, the mixing matrix and monotone modal sources can be precisely obtained. Then, we can estimate the modal parameters by the monotone modal identification method. Sadhu utilized wavelet transform (WT) to transform raw signals into wavelet domain to achieve the sparsity, and then principal component analysis and parallel factor decomposition method were employed for determining the mode shape vectors, natural frequency, and damping, respectively [12, 13]. Yang proposed a scheme using the clustering algorithm to estimate the mode shape matrix in frequency domain and then separated monotone modal responses via linear programming techniques [14]. Yu estimated the mode shape combining with single source point method and time-frequency ratio of mixtures (TIFROM) after executing short-time Fourier transform (STFT), and then natural frequency and damping can be identified by detecting TF ridges [15]. Qin proposed a hybrid method combining [14, 15] to improve the accuracy of matrix estimation and source separation [16]. With more papers being carried out, the SCA method shows more and more powerful ability in modal identification than conventional determined BSS techniques [17, 18].However, the SCA method only works well in the situation of instantaneous mixture, which means that the modal response signals should reach the sensors at the same time [10, 11]. In practical applications, this requirement is also too strict to meet. Therefore, there are two main purposes in this paper: first to provide detailed analysis on the theorem of SCA and point out its limitations in practical applications and then to establish a novel BSS method that can overcome the limitations. With these prime objectives, this paper is organized in the following manner. Section 2 illustrates the detailed theorem of SCA. Section 3 introduces a novel BSS method. Two experimental validations are carried out in Section 4. Then, the conclusion is followed in Section 5.


## 2. Principle of SCA

### 2.1. Motivation of SCA

The SCA framework illustrates that, in the sparse domain (e.g., TF domain or frequency domain), the sparse coefficients of different monotone modals can cluster into the straight line. This phenomenon motivates that it is possible to separate the sparse coefficients based on the clustering features. Therefore, the clustering features in the SCA are the most important issue that should be well explored. Given a numerical example, where three harmonic sourcessit are linearly mixed into two observations xjt:(1)s1t=sin2π1.5t,s2t=sin2π3t,s3t=sin2π5t,(2)x1=s1+0.8s2+0.5s3,x2=0.5s1+0.8s2+s3.We first display the time-series values of two observations by the scatter plot, which is shown in Figure1(a). It can be seen that any useful information on three sources cannot be obtained in time domain. If we use STFT to transform the original signal into TF domain, then we plot the scatter of the real value of two observations, as shown in Figure 1(b). It can be obviously observed that there appear three clustering lines. Meanwhile, the directions of these straight lines correspond to the column vector of the mixing matrix 10.80.50.50.81. By estimating the clustering directions, we can obtain the mixing matrix, and then these sources can be separated by linear programming techniques or L1-norm minimum algorithm. Therefore, it motivates us to first explore the reason that why the clustering features appear.Figure 1
The scatter plot of instantaneous mixture in (a) time domain and (b) TF domain.
(a)
(b)


### 2.2. Lissajous Figure

Before analysis of the SCA, it is necessary to illustrate some essential background knowledge. In mathematics, a Lissajous figure is the graph of a system of parametric equations:(3)x=A⋅sina⋅t,y=B⋅sinb⋅t+δ,which can describe the complex harmonic motion. This motion system was investigated in detail by Lissajous in 1857. The appearance of the figure is highly dependent on the ratio a/b. For a/b=1 and δ=0, the figure is a line, whose direction is determined by A/B, as shown in Figure 2(a). For a/b=1 and δ≠0, the figure is an ellipse, whose shape is determined by A/B and δ, as shown in Figure 2(b).Figure 2
The Lissajous figure: (a) line and (b) ellipse.
(a)
(b)


### 2.3. Sparse Transform Method

Selecting an appropriate sparse method is helpful for executing the SCA algorithm. The time-frequency analysis method has been accepted as the most effective method to achieve sparse in most of papers [10–15, 19–22]. Therefore, we select the STFT as the sparse method in this paper. Given the regular STFT as(4)Gt,ω=∫−∞+∞gu−t⋅su⋅e−iωudu,where gu−t is the moved window and su is the truncated signal. Consider the modified STFT with additional phase shift than regular STFT as(5)Gt,ω=∫−∞+∞gu−t⋅su⋅e−iωu−tdu.According to Parseval’s theorem, the modified STFT can be written as(6)Gt,ω=12π∫−∞+∞g^ω−ξ⋅s^ξ⋅eiξtdξ,where g^ω−ξ is the Fourier transform (FT) of the window function, suppg^⊂−Δ,Δ, and s^ξ is the FT of the signal. Consider a purely harmonic signal st whose frequency is ω0:(7)st=A⋅eiω0t.Considering thats^ξ=2πAδξ−ω0, the STFT of the harmonic signal can be written as(8)Gt,ω=A⋅g^ω−ω0⋅eiω0t.Equation (8) denotes us a first impression on how STFT works [23–27]. The STFT of the harmonic signal is constituted by a series of harmonic signals with the same frequency (which is consistent with the original signal) but different amplitudes (which is determined by both of the signal amplitude and the Fourier transform of window function). Given two harmonic signals with different frequencies:(9)s1t=A1⋅eiω1t,s2t=A2⋅eiω2t.Then, the corresponding STFTs can be written as(10)Gs1t,ω=A1⋅g^ω−ω1⋅eiω1t,Gs2t,ω=A2⋅g^ω−ω2⋅eiω2t.Based on the instantaneous mixture, two sources (9) are first mixed into two observations linearly:(11)x1t=a11⋅s1t+a12⋅s2t,x2t=a21⋅s1t+a22⋅s2t.According to the linearity property of STFT, the STFTs of two observations can be written as(12)Gx1t,ω=a11⋅Gs1t,ω+a12⋅Gs2t,ω=a11⋅A1⋅g^ω−ω1⋅eiω1t+a12⋅A2⋅g^ω−ω2⋅eiω2t,Gx2t,ω=a21⋅Gs1t,ω+a22⋅Gs2t,ω=a21⋅A1⋅g^ω−ω1⋅eiω1t+a22⋅A2⋅g^ω−ω2⋅eiω2t.Herein, it is needed to assumeω1−ω2>2Δ, which means that the frequency distance between two sources is larger than the frequency support of window function. Therefore, this assumption suggests that it is better to select the long window function, which makes sure that the frequency support is narrow enough. Then, we can have the following equations:(13)Gs1t,ω≠0,Gs2t,ω=0,forω∈ω1−Δ,ω1+Δ,(14)Gs1t,ω=0,Gs2t,ω≠0,forω∈ω2−Δ,ω2+Δ.Equations (13) and (14) mean that there exists no overlap between two sources in the TF domain. Therefore, for ω∈ω1−Δ,ω1+Δ, we can construct the following formulation:(15)ReGx1t,ωReGx2t,ω=a11⋅ReGs1t,ωa12⋅ReGs1t,ω=a11⋅A1⋅g^ω−ω1⋅cosω1ta12⋅A1⋅g^ω−ω1⋅cosω1t=a11a12,where Re denotes taking the real part. Then, for ω∈ω2−Δ,ω2+Δ, we can construct the following formulation:(16)ReGx1t,ωReGx2t,ω=a21⋅ReGs2t,ωa22⋅ReGs2t,ω=a21⋅A2⋅g^ω−ω2⋅cosω2ta22⋅A2⋅g^ω−ω2⋅cosω2t=a21a22.Equations (15) and (16) denote that, for the instantaneous mixture, when we plot the scatter of the real part of STFT of two observations, i.e., ReGx1t,ω versus ReGx2t,ω in each frequency bin, it is equal to plot a series of Lissajous figure of the special case (a/b=1, δ=0), i.e., straight lines. These lines have the same direction which is determined by the column vector of mixing matrix, e.g., a11/a12 or a21/a22. This is the reason that why there can appear several clustering lines in the scatter plot (see Figure 1(b)), which can also be employed to estimate the mixing matrix.


### 2.4. Delay Mixture

For the signal of equation (9), we consider the delay mixture:(17)x1t=a11⋅s1t+a12⋅s2t,x2t=a21⋅s1t+∂1+a22⋅s2t+∂2,where ∂i≠0 denotes time delay. According to equation (8), the STFT of two observations can be written as(18)Gx1t,ω=a11⋅Gs1t,ω+a12⋅Gs2t,ω=a11⋅A1⋅g^ω−ω1⋅eiω1t+a12⋅A2⋅g^ω−ω2⋅eiω2t,Gx2t,ω=a21⋅Gs1′t,ω+a22⋅Gs2′t,ω=a21⋅A1⋅g^ω−ω1⋅eiω1t+∂1+a22⋅A2⋅g^ω−ω2⋅eiω2t+∂2.Similar to equation (13), for ω∈ω1−Δ,ω1+Δ, we have(19)Gs1t,ω≠0,Gs1′t,ω≠0,Gs2t,ω=0,Gs2′t,ω=0.Then, we can have the following equation:(20)ReGx1t,ωReGx2t,ω=a11⋅ReGs1t,ωa12⋅ReGs1t,ω=a11⋅cosω1ta12⋅cosω1t+∂1≠a11a12.And forω∈ω2−Δ,ω2+Δ,(21)Gs1t,ω=0,Gs1′t,ω=0,Gs2t,ω≠0,Gs2′t,ω≠0.Then, we can have the following equation:(22)ReGx1t,ωReGx2t,ω=a21⋅ReGs2t,ωa22⋅ReGs2t,ω=a21⋅cosω2ta22⋅cosω2t+∂2≠a21a22.Equations (20) and (22) denote that, for the delay mixture, when we plot the scatter of the real part of STFT of two observations, i.e., ReGx1t,ω versus ReGx2t,ω in each frequency bin, it is equal to plot a series of Lissajous figure with the ellipse case (a/b=1, δ≠0). The shape of these ellipses is determined by both amplitude and time delay in observations. We then utilize the numerical example to further illustrate above analysis. For the numerical example of equation (2), we consider the delay mixture, the scatter plot in the time domain and TF domain is displayed in Figures 3(a) and 3(b), respectively. It can be seen that, both in the time domain and the TF domain, the linear clustering features completely disappear, which lead to that we cannot obtain any useful information on the mixing matrix.Figure 3
The scatter plot of delay mixture in (a) time domain and (b) TF domain.
(a)
(b)From above analysis, it can be known that the sparsity of the observation signals is highly sensitive to time delay. Considering the dynamic engineering in practice, we need to record the response signals by sensors located at different positions. Even if we can measure the signal at the same time, we cannot avoid the propagation delay of sources between all sensors. When the time delay of recorded signals heavily damages the sparsity of the signals, the SCA will fail to separate sources. This drawback of the SCA greatly restricts itself in practical application.


## 2.1. Motivation of SCA

The SCA framework illustrates that, in the sparse domain (e.g., TF domain or frequency domain), the sparse coefficients of different monotone modals can cluster into the straight line. This phenomenon motivates that it is possible to separate the sparse coefficients based on the clustering features. Therefore, the clustering features in the SCA are the most important issue that should be well explored. Given a numerical example, where three harmonic sourcessit are linearly mixed into two observations xjt:(1)s1t=sin2π1.5t,s2t=sin2π3t,s3t=sin2π5t,(2)x1=s1+0.8s2+0.5s3,x2=0.5s1+0.8s2+s3.We first display the time-series values of two observations by the scatter plot, which is shown in Figure1(a). It can be seen that any useful information on three sources cannot be obtained in time domain. If we use STFT to transform the original signal into TF domain, then we plot the scatter of the real value of two observations, as shown in Figure 1(b). It can be obviously observed that there appear three clustering lines. Meanwhile, the directions of these straight lines correspond to the column vector of the mixing matrix 10.80.50.50.81. By estimating the clustering directions, we can obtain the mixing matrix, and then these sources can be separated by linear programming techniques or L1-norm minimum algorithm. Therefore, it motivates us to first explore the reason that why the clustering features appear.Figure 1
The scatter plot of instantaneous mixture in (a) time domain and (b) TF domain.
(a)
(b)


## 2.2. Lissajous Figure

Before analysis of the SCA, it is necessary to illustrate some essential background knowledge. In mathematics, a Lissajous figure is the graph of a system of parametric equations:(3)x=A⋅sina⋅t,y=B⋅sinb⋅t+δ,which can describe the complex harmonic motion. This motion system was investigated in detail by Lissajous in 1857. The appearance of the figure is highly dependent on the ratio a/b. For a/b=1 and δ=0, the figure is a line, whose direction is determined by A/B, as shown in Figure 2(a). For a/b=1 and δ≠0, the figure is an ellipse, whose shape is determined by A/B and δ, as shown in Figure 2(b).Figure 2
The Lissajous figure: (a) line and (b) ellipse.
(a)
(b)


## 2.3. Sparse Transform Method

Selecting an appropriate sparse method is helpful for executing the SCA algorithm. The time-frequency analysis method has been accepted as the most effective method to achieve sparse in most of papers [10–15, 19–22]. Therefore, we select the STFT as the sparse method in this paper. Given the regular STFT as(4)Gt,ω=∫−∞+∞gu−t⋅su⋅e−iωudu,where gu−t is the moved window and su is the truncated signal. Consider the modified STFT with additional phase shift than regular STFT as(5)Gt,ω=∫−∞+∞gu−t⋅su⋅e−iωu−tdu.According to Parseval’s theorem, the modified STFT can be written as(6)Gt,ω=12π∫−∞+∞g^ω−ξ⋅s^ξ⋅eiξtdξ,where g^ω−ξ is the Fourier transform (FT) of the window function, suppg^⊂−Δ,Δ, and s^ξ is the FT of the signal. Consider a purely harmonic signal st whose frequency is ω0:(7)st=A⋅eiω0t.Considering thats^ξ=2πAδξ−ω0, the STFT of the harmonic signal can be written as(8)Gt,ω=A⋅g^ω−ω0⋅eiω0t.Equation (8) denotes us a first impression on how STFT works [23–27]. The STFT of the harmonic signal is constituted by a series of harmonic signals with the same frequency (which is consistent with the original signal) but different amplitudes (which is determined by both of the signal amplitude and the Fourier transform of window function). Given two harmonic signals with different frequencies:(9)s1t=A1⋅eiω1t,s2t=A2⋅eiω2t.Then, the corresponding STFTs can be written as(10)Gs1t,ω=A1⋅g^ω−ω1⋅eiω1t,Gs2t,ω=A2⋅g^ω−ω2⋅eiω2t.Based on the instantaneous mixture, two sources (9) are first mixed into two observations linearly:(11)x1t=a11⋅s1t+a12⋅s2t,x2t=a21⋅s1t+a22⋅s2t.According to the linearity property of STFT, the STFTs of two observations can be written as(12)Gx1t,ω=a11⋅Gs1t,ω+a12⋅Gs2t,ω=a11⋅A1⋅g^ω−ω1⋅eiω1t+a12⋅A2⋅g^ω−ω2⋅eiω2t,Gx2t,ω=a21⋅Gs1t,ω+a22⋅Gs2t,ω=a21⋅A1⋅g^ω−ω1⋅eiω1t+a22⋅A2⋅g^ω−ω2⋅eiω2t.Herein, it is needed to assumeω1−ω2>2Δ, which means that the frequency distance between two sources is larger than the frequency support of window function. Therefore, this assumption suggests that it is better to select the long window function, which makes sure that the frequency support is narrow enough. Then, we can have the following equations:(13)Gs1t,ω≠0,Gs2t,ω=0,forω∈ω1−Δ,ω1+Δ,(14)Gs1t,ω=0,Gs2t,ω≠0,forω∈ω2−Δ,ω2+Δ.Equations (13) and (14) mean that there exists no overlap between two sources in the TF domain. Therefore, for ω∈ω1−Δ,ω1+Δ, we can construct the following formulation:(15)ReGx1t,ωReGx2t,ω=a11⋅ReGs1t,ωa12⋅ReGs1t,ω=a11⋅A1⋅g^ω−ω1⋅cosω1ta12⋅A1⋅g^ω−ω1⋅cosω1t=a11a12,where Re denotes taking the real part. Then, for ω∈ω2−Δ,ω2+Δ, we can construct the following formulation:(16)ReGx1t,ωReGx2t,ω=a21⋅ReGs2t,ωa22⋅ReGs2t,ω=a21⋅A2⋅g^ω−ω2⋅cosω2ta22⋅A2⋅g^ω−ω2⋅cosω2t=a21a22.Equations (15) and (16) denote that, for the instantaneous mixture, when we plot the scatter of the real part of STFT of two observations, i.e., ReGx1t,ω versus ReGx2t,ω in each frequency bin, it is equal to plot a series of Lissajous figure of the special case (a/b=1, δ=0), i.e., straight lines. These lines have the same direction which is determined by the column vector of mixing matrix, e.g., a11/a12 or a21/a22. This is the reason that why there can appear several clustering lines in the scatter plot (see Figure 1(b)), which can also be employed to estimate the mixing matrix.


## 2.4. Delay Mixture

For the signal of equation (9), we consider the delay mixture:(17)x1t=a11⋅s1t+a12⋅s2t,x2t=a21⋅s1t+∂1+a22⋅s2t+∂2,where ∂i≠0 denotes time delay. According to equation (8), the STFT of two observations can be written as(18)Gx1t,ω=a11⋅Gs1t,ω+a12⋅Gs2t,ω=a11⋅A1⋅g^ω−ω1⋅eiω1t+a12⋅A2⋅g^ω−ω2⋅eiω2t,Gx2t,ω=a21⋅Gs1′t,ω+a22⋅Gs2′t,ω=a21⋅A1⋅g^ω−ω1⋅eiω1t+∂1+a22⋅A2⋅g^ω−ω2⋅eiω2t+∂2.Similar to equation (13), for ω∈ω1−Δ,ω1+Δ, we have(19)Gs1t,ω≠0,Gs1′t,ω≠0,Gs2t,ω=0,Gs2′t,ω=0.Then, we can have the following equation:(20)ReGx1t,ωReGx2t,ω=a11⋅ReGs1t,ωa12⋅ReGs1t,ω=a11⋅cosω1ta12⋅cosω1t+∂1≠a11a12.And forω∈ω2−Δ,ω2+Δ,(21)Gs1t,ω=0,Gs1′t,ω=0,Gs2t,ω≠0,Gs2′t,ω≠0.Then, we can have the following equation:(22)ReGx1t,ωReGx2t,ω=a21⋅ReGs2t,ωa22⋅ReGs2t,ω=a21⋅cosω2ta22⋅cosω2t+∂2≠a21a22.Equations (20) and (22) denote that, for the delay mixture, when we plot the scatter of the real part of STFT of two observations, i.e., ReGx1t,ω versus ReGx2t,ω in each frequency bin, it is equal to plot a series of Lissajous figure with the ellipse case (a/b=1, δ≠0). The shape of these ellipses is determined by both amplitude and time delay in observations. We then utilize the numerical example to further illustrate above analysis. For the numerical example of equation (2), we consider the delay mixture, the scatter plot in the time domain and TF domain is displayed in Figures 3(a) and 3(b), respectively. It can be seen that, both in the time domain and the TF domain, the linear clustering features completely disappear, which lead to that we cannot obtain any useful information on the mixing matrix.Figure 3
The scatter plot of delay mixture in (a) time domain and (b) TF domain.
(a)
(b)From above analysis, it can be known that the sparsity of the observation signals is highly sensitive to time delay. Considering the dynamic engineering in practice, we need to record the response signals by sensors located at different positions. Even if we can measure the signal at the same time, we cannot avoid the propagation delay of sources between all sensors. When the time delay of recorded signals heavily damages the sparsity of the signals, the SCA will fail to separate sources. This drawback of the SCA greatly restricts itself in practical application.


## 3. A Novel BSS Method

### 3.1. Theory of the Proposed Method

In this section, we develop an effective BSS method that can overcome the limitations in the SCA method. For a specified harmonic signal whose frequency isω0, we have the reconstruction formulation of STFT as(23)st=ReCg−1∫ω0−Δω0+ΔGt,ωdω,where Cg=1/2∫g^ξdξ. Equation (23) denotes that if we can obtain the frequency band of each source in the TF representation, i.e., ω0−Δ,ω0+Δ, we can reconstruct each source effectively. Therefore, it motivates us to first estimate the frequency band for each source. Considering from the viewpoint of signal energy, at the same time, the energy ratio of the same source recorded by two sensors should be determined, which depends on the source amplitude in recorded signals. Therefore, it inspires us to estimate the frequency band by calculating the TF energy of the signal. Based on the STFT of signals, we first define the formulation of frequency energy function as(24)Eω=∫−∞+∞Gt,ωdt.For the signal of equation (18), we can have the following equation:(25)E1ωE2ω=∫−∞+∞Gx1t,ωdt∫−∞+∞Gx2t,ωdt=∫−∞+∞a11⋅Gs1t,ω+a12⋅Gs2t,ωdt∫−∞+∞a21⋅Gs1′t,ω+a22⋅Gs2′t,ωdt.For the frequency band ofs1t, i.e., ω∈ω1−Δ,ω1+Δ, according to equation (19), equation (25) can be written as(26)E1ωE2ω=∫−∞+∞a11⋅Gs1t,ωdt∫−∞+∞a21⋅Gs1′t,ωdt=∫−∞+∞a11⋅A1⋅g^ω−ω1⋅eiω1tdt∫−∞+∞a21⋅A1⋅g^ω−ω1⋅eiω1t+∂1dt=a11a21⋅A1⋅g^ω−ω1A1⋅g^ω−ω1⋅∫−∞+∞eiω1tdt∫−∞+∞eiω1t+∂1dt=a11a21,where eiω1t=1. Similarly, for the frequency band of s2t, i.e., ω∈ω2−Δ,ω2+Δ, equation (25) can be written as(27)E1ωE2ω=∫−∞+∞a12⋅Gs2t,ωdt∫−∞+∞a22⋅Gs2′t,ωdt=∫−∞+∞a12⋅A2⋅g^ω−ω2⋅eiω2tdt∫−∞+∞a22⋅A2⋅g^ω−ω2⋅eiω2t+∂2dt=a12a22⋅A2⋅g^ω−ω2A2⋅g^ω−ω2⋅∫−∞+∞eiω2tdt∫−∞+∞eiω2t+∂2dt=a12a22.Equations (26) and (27) denote that, for the frequency band of each source, the ratio of frequency energy is consistent, which is determined by the source amplitude in recorded signals. When we plot the scatter of the frequency energy data of two observations, i.e. E1ω versus E2ω, if a11/a21≠a12/a22, there will appear several clustering lines which can correspond to the frequency band of each source. Therefore, we can estimate the frequency band of each source according to these clustering features. For the delay mixture of numerical signals of equation (2), the scatter plot of frequency energy is shown as Figure 4(a). It can be seen that there are three obvious clustering lines.Figure 4
(a) The scatter plot of frequency energy data; (b) the sum of frequency energy data.
(a)
(b)In Figure4(a), each point in the clustering lines corresponds to the frequency bin of source in TF representation. Then, it is required to develop an effective procedure to estimate these points. As observing the scatter plot, each clustering line has the maximum value located at the end of line (as red label “∗”), and other points have the same cosine distance with this end point. The point of maximum value can be detected by the peak data by summing all frequency energy data:(28)Eω=∑i=1mEiω.Therefore, we can first obtain the maximum value point, as shown in Figure4(b) (as red label “∗”), and then calculate the cosine distance with it to estimate other points by(29)CDE1ω,E2ω,…,Enω,E1ωp,E2ωp,…,Enωp<ε,where ωp is the detected peak frequency and ε is a low value which is suggested to be empirically set as 0.004. If condition (29) is satisfied, these scatter points should belong to the same source. By repeating this procedure, we can obtain the frequency bands of all sources. After frequency band being estimated, we can reconstruct all sources by inversion of STFT. For the numerical signal, the three recovered sources are shown in Figure 5. It can be observed that three sources are well separated as monocomponents.Figure 5
The recovered three sources S1, S2, and S3.A known drawback of STFT reconstruction is the boundary effect, i.e., the amplitude in the beginning and end of the sources cannot be well recovered. In order to eliminate the boundary effect of the recovered sources, we introduce a padding line, as(30)Lt=∫−∞+∞gudu∫t1t2gu−tdu,where gu−t is the moved window and t1,t2 denotes the discrete beginning and end time point of the analysed signal. The padding line of the numerical signal is shown in Figure 6. And then, each source is corrected through multiplying this padding line. The three sources that consider padding method are shown in Figure 7. It can be seen that the amplitudes of all sources are also recovered precisely.Figure 6
The padding line.Figure 7
The recovered three sources considering padding line S1, S2, and S3.Therefore, the proposed BSS method can be summarized as follows:(1)
Utilize STFT to transform recorded signals into the TF domain(2)
Obtain the frequency energy data according to equation (24)(3)
Sum all frequency energy data and pick the peak data by the peak detection method(4)
Calculate the cosine distance of other points with detected peak data to estimate the frequency band(5)
Recover each source according to the estimated frequency band by inversion of STFT(6)
Eliminate the boundary effect of each source by the padding method


### 3.2. Numerical Signal Analysis

In this subsection, we consider the performance of the proposed method in analysing a five-component signal added with heavy Gaussian noise. This five-component signal mixed into two observations is modelled as(31)S1t=sin2π30t+sin2π70t+sin2π110t+sin2π150t+sin2π190t,S2t=0.5sin2π30t+0.9sin2π70t+1.3sin2π110t+1.7sin2π150t+2.2sin2π190t,where the sampling frequency is 512 Hz and sampling time is 2 s. The waveforms of the two observations are plotted in Figure 8. In Figure 9, the time-frequency representations of two observations generated by STFT are displayed. Then, the scatter of the frequency energy data is calculated and plotted in Figure 10. It can be seen that there are five obvious clustering lines in Figure 10, which correspond to the frequency band of five sources. The eventual separated sources are plotted in Figure 11. It is shown that all sources are well separated into monocomponent signals.Figure 8
Two observations: (a) S1 and (b) S2.
(a)
(b)Figure 9
Time-frequency representation of (a) S1 and (b) S2.
(a)
(b)Figure 10
The scatter plot of frequency energy data.Figure 11
Five separated sources, from low-frequency components to high-frequency components (from top to bottom).


## 3.1. Theory of the Proposed Method

In this section, we develop an effective BSS method that can overcome the limitations in the SCA method. For a specified harmonic signal whose frequency isω0, we have the reconstruction formulation of STFT as(23)st=ReCg−1∫ω0−Δω0+ΔGt,ωdω,where Cg=1/2∫g^ξdξ. Equation (23) denotes that if we can obtain the frequency band of each source in the TF representation, i.e., ω0−Δ,ω0+Δ, we can reconstruct each source effectively. Therefore, it motivates us to first estimate the frequency band for each source. Considering from the viewpoint of signal energy, at the same time, the energy ratio of the same source recorded by two sensors should be determined, which depends on the source amplitude in recorded signals. Therefore, it inspires us to estimate the frequency band by calculating the TF energy of the signal. Based on the STFT of signals, we first define the formulation of frequency energy function as(24)Eω=∫−∞+∞Gt,ωdt.For the signal of equation (18), we can have the following equation:(25)E1ωE2ω=∫−∞+∞Gx1t,ωdt∫−∞+∞Gx2t,ωdt=∫−∞+∞a11⋅Gs1t,ω+a12⋅Gs2t,ωdt∫−∞+∞a21⋅Gs1′t,ω+a22⋅Gs2′t,ωdt.For the frequency band ofs1t, i.e., ω∈ω1−Δ,ω1+Δ, according to equation (19), equation (25) can be written as(26)E1ωE2ω=∫−∞+∞a11⋅Gs1t,ωdt∫−∞+∞a21⋅Gs1′t,ωdt=∫−∞+∞a11⋅A1⋅g^ω−ω1⋅eiω1tdt∫−∞+∞a21⋅A1⋅g^ω−ω1⋅eiω1t+∂1dt=a11a21⋅A1⋅g^ω−ω1A1⋅g^ω−ω1⋅∫−∞+∞eiω1tdt∫−∞+∞eiω1t+∂1dt=a11a21,where eiω1t=1. Similarly, for the frequency band of s2t, i.e., ω∈ω2−Δ,ω2+Δ, equation (25) can be written as(27)E1ωE2ω=∫−∞+∞a12⋅Gs2t,ωdt∫−∞+∞a22⋅Gs2′t,ωdt=∫−∞+∞a12⋅A2⋅g^ω−ω2⋅eiω2tdt∫−∞+∞a22⋅A2⋅g^ω−ω2⋅eiω2t+∂2dt=a12a22⋅A2⋅g^ω−ω2A2⋅g^ω−ω2⋅∫−∞+∞eiω2tdt∫−∞+∞eiω2t+∂2dt=a12a22.Equations (26) and (27) denote that, for the frequency band of each source, the ratio of frequency energy is consistent, which is determined by the source amplitude in recorded signals. When we plot the scatter of the frequency energy data of two observations, i.e. E1ω versus E2ω, if a11/a21≠a12/a22, there will appear several clustering lines which can correspond to the frequency band of each source. Therefore, we can estimate the frequency band of each source according to these clustering features. For the delay mixture of numerical signals of equation (2), the scatter plot of frequency energy is shown as Figure 4(a). It can be seen that there are three obvious clustering lines.Figure 4
(a) The scatter plot of frequency energy data; (b) the sum of frequency energy data.
(a)
(b)In Figure4(a), each point in the clustering lines corresponds to the frequency bin of source in TF representation. Then, it is required to develop an effective procedure to estimate these points. As observing the scatter plot, each clustering line has the maximum value located at the end of line (as red label “∗”), and other points have the same cosine distance with this end point. The point of maximum value can be detected by the peak data by summing all frequency energy data:(28)Eω=∑i=1mEiω.Therefore, we can first obtain the maximum value point, as shown in Figure4(b) (as red label “∗”), and then calculate the cosine distance with it to estimate other points by(29)CDE1ω,E2ω,…,Enω,E1ωp,E2ωp,…,Enωp<ε,where ωp is the detected peak frequency and ε is a low value which is suggested to be empirically set as 0.004. If condition (29) is satisfied, these scatter points should belong to the same source. By repeating this procedure, we can obtain the frequency bands of all sources. After frequency band being estimated, we can reconstruct all sources by inversion of STFT. For the numerical signal, the three recovered sources are shown in Figure 5. It can be observed that three sources are well separated as monocomponents.Figure 5
The recovered three sources S1, S2, and S3.A known drawback of STFT reconstruction is the boundary effect, i.e., the amplitude in the beginning and end of the sources cannot be well recovered. In order to eliminate the boundary effect of the recovered sources, we introduce a padding line, as(30)Lt=∫−∞+∞gudu∫t1t2gu−tdu,where gu−t is the moved window and t1,t2 denotes the discrete beginning and end time point of the analysed signal. The padding line of the numerical signal is shown in Figure 6. And then, each source is corrected through multiplying this padding line. The three sources that consider padding method are shown in Figure 7. It can be seen that the amplitudes of all sources are also recovered precisely.Figure 6
The padding line.Figure 7
The recovered three sources considering padding line S1, S2, and S3.Therefore, the proposed BSS method can be summarized as follows:(1)
Utilize STFT to transform recorded signals into the TF domain(2)
Obtain the frequency energy data according to equation (24)(3)
Sum all frequency energy data and pick the peak data by the peak detection method(4)
Calculate the cosine distance of other points with detected peak data to estimate the frequency band(5)
Recover each source according to the estimated frequency band by inversion of STFT(6)
Eliminate the boundary effect of each source by the padding method


## 3.2. Numerical Signal Analysis

In this subsection, we consider the performance of the proposed method in analysing a five-component signal added with heavy Gaussian noise. This five-component signal mixed into two observations is modelled as(31)S1t=sin2π30t+sin2π70t+sin2π110t+sin2π150t+sin2π190t,S2t=0.5sin2π30t+0.9sin2π70t+1.3sin2π110t+1.7sin2π150t+2.2sin2π190t,where the sampling frequency is 512 Hz and sampling time is 2 s. The waveforms of the two observations are plotted in Figure 8. In Figure 9, the time-frequency representations of two observations generated by STFT are displayed. Then, the scatter of the frequency energy data is calculated and plotted in Figure 10. It can be seen that there are five obvious clustering lines in Figure 10, which correspond to the frequency band of five sources. The eventual separated sources are plotted in Figure 11. It is shown that all sources are well separated into monocomponent signals.Figure 8
Two observations: (a) S1 and (b) S2.
(a)
(b)Figure 9
Time-frequency representation of (a) S1 and (b) S2.
(a)
(b)Figure 10
The scatter plot of frequency energy data.Figure 11
Five separated sources, from low-frequency components to high-frequency components (from top to bottom).


## 4. Experimental Validation

In this section, we employ two experiments to validate the proposed method. These two experiments have been analysed in two published papers [15, 17]. Therefore, their analysed results can be the comparison references. To test the ability of the proposed method comparable to the SCA, we also display the results analysed by SCA.


### 4.1. Experiment 1

As shown in Figure12, the structure is a uniform TC4 titanium-alloy column of 0.38∗0.038∗0.006m3, and three acceleration sensors are mounted on the beam which is excited by impact hammer. The sampling frequency is 2560 Hz, and we take 0.5 s data to analyse. The acceleration responses X1, X2, and X3 are shown in Figure 13. For comparative analysis, we list the identified results in [15], such as frequency and damping ratio, as shown in Table 1.Figure 12
The experimental structure in [15].Figure 13
The time waveform and the spectrum of the recorded vibration signals X1, X2, and X3.Table 1
Identified results of the structure in experiment 1.
Modal
Reference [15]
SCA
Proposed method
Frequency (Hz)
1st
31.28
30.72
30.78
2nd
216.5
215.6
215.6
3rd
596.8
585.9
585.9
4th
1115
1115.2
1115.3
Dam (%)
1st
2.661
2.733
2.863
2nd
0.511
0.508
0.505
3rd
0.601
0.606
0.585
4th
0.242
0.234
0.234To show the processing procedure of the SCA method, we first utilize STFT to transform three recorded signals into the TF domain. Then, the scatter of original time-series signals and the real part of the STFT results are plotted in Figures14(a) and 14(b), respectively. It can be seen that the original signal in time domain cannot provide any useful information about the mixing matrix. However, in Figure 14(b), there obviously appear four clustering lines for the recorded signals, which should correspond to four column vectors of the mixing matrix. Then, we can separate each source using this mixing matrix, as shown in Figure 15. It can be seen that, through the SCA method, each source is well separated as monotone exponentially decaying sinusoid. Then, we utilize the monotone mode method to identify modal parameters, and the identification results are listed in Table 1. It is known that the SCA provides a similar modal identification result with [15].Figure 14
The scatter plot in (a) time domain and (b) TF domain generated by SCA.
(a)
(b)Figure 15
The sources S1, S2, S3, and S4 and their spectrums separated by SCA.Then, for the proposed method, according to equation (24), we first calculate the frequency energy function. The detected peaks of the sum of Eiω are shown in Figure 16(a), where “∗” denotes the detected peaks. It can be seen that six peaks are detected. Among them, two more sources are discovered than SCA, which is caused by unexpected interference. In Figure 16(b), we plot the scatter of frequency energy data, which also shows six obvious clustering lines. To facilitate comparison with SCA, we do not list the separated results of two interference sources. Therefore, four separated sources are shown in Figure 17. It is shown that each source is well recovered as the monomodal response. The estimated frequencies and damping are listed in Table 1, which show satisfactory accuracy [15] as well. It can be concluded that when the sparsity of signal in the TF domain is satisfactory, the monomodal sources can be well separated both by SCA and the proposed method.Figure 16
(a) The sum of frequency energy data and (b) the scatter plot of the frequency energy data generated by the proposed method.
(a)
(b)Figure 17
The sources S1, S2, S3, and S4 and their spectrums separated by the proposed method.


### 4.2. Experiment 2

As shown in Figure18, the structure is a uniform steel cantilever beam of 0.9∗0.05∗0.008m3, and five displacement sensors are settled upon the beam which is excited by impact hammer. The sampling frequency is 1600 Hz, and we take 1000 samples to analyse. The vibration responses X1–X5 are shown in Figure 19. Meanwhile, we first list the identified modal parameters by the method in reference [17] in Table 2.Figure 18
The experimental structure in [17].Figure 19
The time waveform and the spectrum of the recorded vibration signals X1, X2, X3, X4, and X5.Table 2
Identified results of the structure in experiment 2.
Mode
Reference [17]
SCA
Proposed method
Frequency (Hz)
1st
10.07
8.948
8.957
2nd
56.27
55.87
55.85
3rd
155.7
156.5
156.4
4th
304.2
306.1
306.1
5th
500.9
481.4
505.4To illustrate the ability of the proposed method applied in the underdetermined case, we only deal with three sensor data, i.e., X1, X2, and X5. The scatter plot of STFT results of the recorded signals is first displayed in Figure20. It can be seen that, due to influence of time delay, there appear several ellipses, which make the SCA hard to estimate the accurate mixing matrix. The separated results by the SCA are also displayed in Figure 21. It can be observed that the recovered sources S3, S4, and S5 are unwillingly mixed with other sources.Figure 20
The scatter plot generated by the SCA.Figure 21
The separated sources S1, S2, S3, S4, and S5 and their spectrums by the SCA.For our proposed method, according to equation (24), we first calculate the frequency energy Eiω. And the scatter plot of frequency energy data is shown in Figure 22, where “∗” denotes the detected peaks. It can be seen that there appear five obvious clustering lines, which correspond to the frequency band of five sources. The sources separated by the proposed method are listed in Figure 23, which show all monomodal responses being recovered successfully. The identified parameters are listed in Table 2. Due to the recovered sources being close to harmonic signals, it cannot provide accurate estimation for damping, which is neglected in the identified parameters. From the identified parameters, it can be seen that the identified modal parameters show satisfactory accuracy with [17].Figure 22
The scatter plot of the frequency energy data generated by the proposed method.Figure 23
The sources separated S1, S2, S3, S4, and S5 and their spectrums by the proposed method.From above analysis, it can be known that, when the sparsity of the signal in the TF domain is heavily damaged by time delay, the SCA may fail to achieve accurate estimation of mixing matrix and recovery of sources. However, our proposed method can deal with it successfully. Furthermore, we consider the computational cost of the proposed method in dealing with the experimental data. Because it is helpful for judging whether the proposed method can be used in real-time applications or not. The tested computer configuration is as follows: Intel Core i7-6500 2.5 GHz, 8.0 GB of DDR3 RAM, Windows 10 OS, and MATLAB version R2016a. The required computation times in the first experiment and the second experiment are 0.91 s and 0.82 s, respectively. It can be concluded that the proposed methods can finish the processing within one second, which is in the acceptable range in real-time applications.


## 4.1. Experiment 1

As shown in Figure12, the structure is a uniform TC4 titanium-alloy column of 0.38∗0.038∗0.006m3, and three acceleration sensors are mounted on the beam which is excited by impact hammer. The sampling frequency is 2560 Hz, and we take 0.5 s data to analyse. The acceleration responses X1, X2, and X3 are shown in Figure 13. For comparative analysis, we list the identified results in [15], such as frequency and damping ratio, as shown in Table 1.Figure 12
The experimental structure in [15].Figure 13
The time waveform and the spectrum of the recorded vibration signals X1, X2, and X3.Table 1
Identified results of the structure in experiment 1.
Modal
Reference [15]
SCA
Proposed method
Frequency (Hz)
1st
31.28
30.72
30.78
2nd
216.5
215.6
215.6
3rd
596.8
585.9
585.9
4th
1115
1115.2
1115.3
Dam (%)
1st
2.661
2.733
2.863
2nd
0.511
0.508
0.505
3rd
0.601
0.606
0.585
4th
0.242
0.234
0.234To show the processing procedure of the SCA method, we first utilize STFT to transform three recorded signals into the TF domain. Then, the scatter of original time-series signals and the real part of the STFT results are plotted in Figures14(a) and 14(b), respectively. It can be seen that the original signal in time domain cannot provide any useful information about the mixing matrix. However, in Figure 14(b), there obviously appear four clustering lines for the recorded signals, which should correspond to four column vectors of the mixing matrix. Then, we can separate each source using this mixing matrix, as shown in Figure 15. It can be seen that, through the SCA method, each source is well separated as monotone exponentially decaying sinusoid. Then, we utilize the monotone mode method to identify modal parameters, and the identification results are listed in Table 1. It is known that the SCA provides a similar modal identification result with [15].Figure 14
The scatter plot in (a) time domain and (b) TF domain generated by SCA.
(a)
(b)Figure 15
The sources S1, S2, S3, and S4 and their spectrums separated by SCA.Then, for the proposed method, according to equation (24), we first calculate the frequency energy function. The detected peaks of the sum of Eiω are shown in Figure 16(a), where “∗” denotes the detected peaks. It can be seen that six peaks are detected. Among them, two more sources are discovered than SCA, which is caused by unexpected interference. In Figure 16(b), we plot the scatter of frequency energy data, which also shows six obvious clustering lines. To facilitate comparison with SCA, we do not list the separated results of two interference sources. Therefore, four separated sources are shown in Figure 17. It is shown that each source is well recovered as the monomodal response. The estimated frequencies and damping are listed in Table 1, which show satisfactory accuracy [15] as well. It can be concluded that when the sparsity of signal in the TF domain is satisfactory, the monomodal sources can be well separated both by SCA and the proposed method.Figure 16
(a) The sum of frequency energy data and (b) the scatter plot of the frequency energy data generated by the proposed method.
(a)
(b)Figure 17
The sources S1, S2, S3, and S4 and their spectrums separated by the proposed method.


## 4.2. Experiment 2

As shown in Figure18, the structure is a uniform steel cantilever beam of 0.9∗0.05∗0.008m3, and five displacement sensors are settled upon the beam which is excited by impact hammer. The sampling frequency is 1600 Hz, and we take 1000 samples to analyse. The vibration responses X1–X5 are shown in Figure 19. Meanwhile, we first list the identified modal parameters by the method in reference [17] in Table 2.Figure 18
The experimental structure in [17].Figure 19
The time waveform and the spectrum of the recorded vibration signals X1, X2, X3, X4, and X5.Table 2
Identified results of the structure in experiment 2.
Mode
Reference [17]
SCA
Proposed method
Frequency (Hz)
1st
10.07
8.948
8.957
2nd
56.27
55.87
55.85
3rd
155.7
156.5
156.4
4th
304.2
306.1
306.1
5th
500.9
481.4
505.4To illustrate the ability of the proposed method applied in the underdetermined case, we only deal with three sensor data, i.e., X1, X2, and X5. The scatter plot of STFT results of the recorded signals is first displayed in Figure20. It can be seen that, due to influence of time delay, there appear several ellipses, which make the SCA hard to estimate the accurate mixing matrix. The separated results by the SCA are also displayed in Figure 21. It can be observed that the recovered sources S3, S4, and S5 are unwillingly mixed with other sources.Figure 20
The scatter plot generated by the SCA.Figure 21
The separated sources S1, S2, S3, S4, and S5 and their spectrums by the SCA.For our proposed method, according to equation (24), we first calculate the frequency energy Eiω. And the scatter plot of frequency energy data is shown in Figure 22, where “∗” denotes the detected peaks. It can be seen that there appear five obvious clustering lines, which correspond to the frequency band of five sources. The sources separated by the proposed method are listed in Figure 23, which show all monomodal responses being recovered successfully. The identified parameters are listed in Table 2. Due to the recovered sources being close to harmonic signals, it cannot provide accurate estimation for damping, which is neglected in the identified parameters. From the identified parameters, it can be seen that the identified modal parameters show satisfactory accuracy with [17].Figure 22
The scatter plot of the frequency energy data generated by the proposed method.Figure 23
The sources separated S1, S2, S3, S4, and S5 and their spectrums by the proposed method.From above analysis, it can be known that, when the sparsity of the signal in the TF domain is heavily damaged by time delay, the SCA may fail to achieve accurate estimation of mixing matrix and recovery of sources. However, our proposed method can deal with it successfully. Furthermore, we consider the computational cost of the proposed method in dealing with the experimental data. Because it is helpful for judging whether the proposed method can be used in real-time applications or not. The tested computer configuration is as follows: Intel Core i7-6500 2.5 GHz, 8.0 GB of DDR3 RAM, Windows 10 OS, and MATLAB version R2016a. The required computation times in the first experiment and the second experiment are 0.91 s and 0.82 s, respectively. It can be concluded that the proposed methods can finish the processing within one second, which is in the acceptable range in real-time applications.


## 5. Conclusions

In this paper, we focus on a hot topic on application of underdetermined BSS in structural dynamic analysis. The SCA technique, as an effective underdetermined BSS method, has drawn many attentions in recent papers. However, all of these papers only pay attention to the application of SCA, but none of them illustrates its theorem in detail, which lead to some limitations of SCA being neglected. Therefore, we start with the classical TF method and harmonic signal, to illustrate the detailed theorem of SCA. Meanwhile, we point out its limitation in practical engineering. Furthermore, we propose a novel BSS method. By means of TF transform and estimation of frequency band, the proposed method can effectively deal with the problem existed in SCA.The time delay between recorded signals is mainly determined by the propagation speed and distance. For the structures in experiments 1 and 2, they have the similar width and height, but distinct lengths. The structure in experiment 2 is obviously longer than the structure in experiment 1. Under the similar condition of propagation speed, the structure length has become the critical influence of the time delay. It is obvious that a longer structure will produce a larger time delay. Therefore, we suggest that the SCA technique should be applied to analyse small structures, which makes sure the time delay being in the acceptable range.



---

*Source: 1637163-2019-10-08.xml*