A Novel Underdetermined Blind Source Separation Method and Its Application to Source Contribution Quantitative Estimation

Abstract

To identify the major vibration and radiation noise, a source contribution quantitative estimation method is proposed based on underdetermined blind source separation. First, the single source points (SSPs) are identified by directly searching the identical normalized time-frequency vectors of mixed signals, which can improve the efficiency and accuracy in identifying SSPs. Then, the mixing matrix is obtained by hierarchical clustering, and source signals can also be recovered by the least square method. Second, the optimal combination coefficients between source signals and mixed signals can be calculated based on minimum redundant error energy. Therefore, mixed signals can be optimally linearly combined by source signals via the coefficients. Third, the energy elimination method is used to quantitatively estimate source contributions. Finally, the effectiveness of the proposed method is verified via numerical case studies and experiments with a cylindrical structure, and the results show that source signals can be effectively recovered, and source contributions can be quantitatively estimated by the proposed method.

1. Introduction

Vibration and radiation noise have a significant effect on the safety and stability of some mechanical systems [1,2], for example, excessive noise of underwater vehicle will interfere with its own detection accuracy. Independently acquiring information from each source of the mechanical system can help to quickly judge its running state. However, in practice, the information measured by sensors is the superposition of some sources, because different components of the mechanical system will interfere with each other, which makes it difficult to directly measure the source information [3]. Therefore, some supplementary signal processing methods are needed to further process the collected information to obtain the expected source signals [1]. Among post-processing approaches, blind source separation (BSS) has demonstrated its usefulness in separating sources from mixed signals, due to its simplicity and effectiveness. More importantly, BSS can be utilized without the structure models and the transmission paths that are difficult to be obtained, and therefore BSS has been widely used in practice [4,5,6,7,8,9]. However, most of these methods are mainly designed for (over)determined BSS where the number of sensors is no smaller than that of sources, and thus they may fail when dealing with underdetermined cases. Therefore, we mainly address the problem of underdetermined BSS (UBSS) where the number of mixed signals is smaller than that of sources in this study. In addition, reducing the vibration of the major sources rather than all sources can achieve satisfactory results with smaller cost [1,10]. Therefore, how to evaluate source contribution quantitatively is another problem to be addressed in our study.

Recently, UBSS has gained considerable research interests and various excellent methods have been proposed to tackle this problem. The sparsity of source signals is often exploited for UBSS. Bofill and Zibulevsky proposed a two-step algorithm, i.e., mixing matrix estimation and source recovery, that exploits the sparsity of source signals in the frequency domain [11]. The framework of the two-step algorithm is used by many researchers later. Wigner-Ville distribution (WVD) is usually adopted to transform mixed signals from the time domain to time-frequency (TF) domain, which can greatly increase the sparsity of source signals. Aïssa-El-Bey et al. proposed an effective WVD-based UBSS method that does not require disjoint source signals in the whole TF domain and only needs some TF points where only one source exists, called single source points (SSPs) [12]. After obtaining the mixing matrix, source signals are recovered by the subspace method under the assumption that the number of active source signals must be smaller than sensors at any TF point. Then, this restriction is relaxed by Peng and Xiang, which only requires that the number of active source signals at any TF point should not exceed that of sensors [13]. Xie et al. further relaxed the sparsity constraint and developed a new WVD-based UBSS algorithm to recover the sources exactly at every auto-term TF point, provided that the number of the sources is smaller than the twice of the sensors [14]. Though WVD can provide a high concentration of signal energy in the TF domain, it has some inherent disadvantages of cross-term disturbance and high computational complexity. Therefore, short-time Fourier transform (STFT) is also adopted by many researchers recently [15,16], due to its freedom of cross terms and its easiness to be inverted. Reju et al. proposed a mixing matrix estimation algorithm based on STFT and proved that a TF point is an SSP when the absolute direction of the real and the imaginary parts of TF vectors of the mixed signals are the same [17]. Li et al. also proposed an effective algorithm to identify SSPs by utilizing time-frequency coefficients of the mixed signals and the complex conjugates of the coefficients [18]. However, most of the SSPs-based UBSS methods focus merely on the nature of single SSP ignoring the relationship between SSPs, which leads to low identification accuracy of SSPs, especially in noisy cases. Zhen et al. proposed a UBSS method based on STFT in which SSPs are identified by sparse coding. Since the sparse coding strategy considers the linear relations between SSPs, Zhen’s method can obtain excellent estimation performance even in low signal to noise ratio (SNR) cases [19]. However, the SSPs identification effectiveness of Zhen’s method is still low when SSPs are rare. Furthermore, it is a great challenge to select the regularization parameter, since improper parameter will lead to larger errors of the estimated sources. Hence, it is still difficult to identify SSPs accurately and efficiently.

After source signals are recovered, source contributions can be quantitatively estimated. According to the percentage of the source contribution, the major sources to the measurement points can also be identified. Cheng et al. proposed a source contribution evaluation method of mechanical vibration signals via enhanced independent component analysis for complete BSS [4,11]. Zhang et al. proposed a source contribution method with kurtosis-based constrained independent component analysis, which can quantitatively estimate source contribution by reducing the energy of the extracted independent components from mixed signals in each extraction [1]. However, the separation matrix should be an orthogonal matrix, which cannot be satisfied in UBSS. In this paper, we are mainly addressing the issue of source contribution estimation for UBSS.

To accurately recover source signals and identify the major sources, an effective UBSS method and a source contribution quantitative estimation method are proposed, in which SSPs can be identified by directly searching the identical normalized TF vectors of mixed signals instead of the sparse coding in Zhen’s method. After source signals are well recovered by the proposed UBSS method, the optimal combination coefficients between source signals and mixed signals can be obtained based on minimum redundant error energy. Then, the contributions of the recovered source signals are quantitatively estimated by using energy elimination method. Finally, the effectiveness of the proposed method is verified by some numerical simulations and experimental studies on a cylindrical shell structure.

The organization of this paper is as follows. In Section 2, an effective UBSS method is proposed. In Section 3, the source contribution quantitative estimation method is given. In Section 4, the proposed UBSS method and the contribution quantitative estimation method are verified by the numerical case studies. In Section 5, the proposed method is further verified with experimental studies on a cylindrical shell structure. Some conclusions are drawn in Section 6.

2. Underdetermined Blind Source Separation

2.1. Basic Theory

The linear instantaneous mixed model of UBSS can be expressed as

$$\mathbf{x}(t)=\mathbf{A}\mathbf{s}(t),$$

(1)

where $\mathbf{x}(t)={[x_1(t),\hspace{0.17em}{x}_2(t),\hspace{0.17em}\dots ,x_N(t)]}^{\mathrm{T}}$ and $\mathbf{s}(t)={[s_1(t),\hspace{0.17em}{s}_2(t),\hspace{0.17em}\dots ,s_M(t)]}^{\mathrm{T}}$ are the mixed vector and the source vector in the time domain, respectively, and ${·}^{\mathrm{T}}$ represents the transpose operation; N and M (N < M) are the number of mixed signals and source signals, respectively; $\mathbf{A}=[{\mathbf{a}}_1,\hspace{0.17em}{\mathbf{a}}_2,\hspace{0.17em}\dots ,\hspace{0.17em}{\mathbf{a}}_M]$ is the mixing matrix with ${\mathbf{a}}_i$ as its column. The aim of UBSS is to estimate source signals without any prior information of $s(t)$ or $\mathbf{A}$, except that N < M.

To increase the sparsity of source signals, the above linear instantaneous mixed model can be transformed into the TF domain as Equation (2) or Equation (3) by STFT.

$$\mathbf{X}\left(t,f\right)=AS\left(t,f\right),$$

(2)

$$\mathbf{X}\left(t,f\right)={\displaystyle \sum _{i=1}^M{a}_i{S}_i\left(t,f\right)},$$

(3)

where $X\left(t,f\right)={\left[X_1\left(t,f\right),\hspace{0.17em}{X}_2\left(t,f\right),\hspace{0.17em}\dots ,X_N\left(t,f\right)\right]}^{\mathrm{T}}$ and $S\left(t,f\right)={\left[S_1\left(t,f\right),\hspace{0.17em}{S}_2\left(t,f\right),\dots ,S_M\left(t,f\right)\right]}^{\mathrm{T}}$ are the STFT coefficients of $\mathbf{x}\left(t\right)$ and $\mathbf{s}\left(t\right)$ at TF point $\left(t,f\right)$, respectively.

2.2. Proposed Mixing Matrix Estimation Method

The ideal goal of UBSS is to estimate source signals without any prior information, except that $N<M$. Actually, it seems almost impossible to obtain an effective estimation of source signals if we know nothing about $s(t)$ or $\mathbf{A}$. Therefore, some assumptions are given first.

Assumption 1.

For each source signal $s_i(t)$, there are some TF points $(t,f)$ where only $S_i(t,f)$ is dominant, i.e., $\left|S_i(t,f)\right|\gg \left|S_j(t,f)\right|$, $\forall \hspace{0.17em}j\ne i$.

Assumption 2.

Source signals are mutually independent.

Assumption 3.

Any $N\times N$ submatrix of the mixing matrix $A$ is of full rank.

These three assumptions could hold in many practical cases and have been widely used in recent UBSS methods. Assumption 1 is used to guarantee the existence of SSPs in the process of mixing matrix estimation. However, Assumption 1 is not necessary for recovering source signals, that is, if the mixing matrix is known or can be estimated by other methods, Assumption 1 can be removed. Assumption 2 and Assumption 3 are used to increase the stability of the SSPs identification method. Assumption 2 is also used to guarantee that there is no cross energy among source signals, which will be used in the source contribution estimation. Besides, Assumption 3 is also used to guarantee that all source signals could be correctly recovered.

Now, at any TF point, say $(u,v)$, if only one source is active, say $S_i(u,v)$, i.e., $(u,v)$ is an SSP corresponding to $s_i(t)$, then Equation (3) can be rewritten as

$$\mathbf{X}(u,v)=S_i(u,v)a_i.$$

(4)

Equation (4) shows that the TF vector of mixed signals at TF point $\left(u,v\right)$ is collinear with the $i$-th column of the mixing matrix. It can be also obtained that all TF vectors of mixed signals at SSPs corresponding to $s_i(t)$ will be collinear with $a_i$, that is, all SSPs corresponding to the same active source could be linearly represented by each other. Assume that $S_i\left(\psi ,\omega \right)$ is also an SSP corresponding to $s_i(t)$, then we will obtain

$$\mathbf{X}(u,v)=r\mathbf{X}(\psi ,\omega ),$$

(5)

where $r$ is a real coefficient.

Now, the problem is how to identify the TF vectors that satisfy Equation (5) from all TF vectors of mixed signals.

From Equation (4), $\mathbf{X}\left(u,v\right)$ can be normalized as

$$\begin{array}{ll}\tilde{\mathbf{X}}\left(u,v\right)& =\frac{S_i\left(u,v\right)a_i}{{\Vert S_i\left(u,v\right)a_i\Vert }_2}\\ & =\frac{a_i}{{\Vert a_i\Vert }_2}\\ & ={\tilde{a}}_i\end{array},$$

(6)

where ῀ represents the normalized vector and ${\Vert \Vert }_2$ is the 2-norm. Similarly, the normalized vector of $\mathbf{X}\left(\psi ,\omega \right)$ can also be written as

$$\tilde{\mathbf{X}}\left(\psi ,\omega \right)={\tilde{a}}_i$$

(7)

As shown in Equations (6) and (7), all the normalized TF vectors at the SSPs corresponding to the i-th dominant source will be equal to the normalized vector of $a_i$. Therefore, SSPs can be identified by searching the identical normalized TF vectors of mixed signals, i.e.,

$$\tilde{\mathbf{X}}(u,v)=\tilde{\mathbf{X}}(\psi ,\omega )$$

(8)

Equations (5) and (8) are equivalent. Therefore, SSPs can be identified by checking whether normalized TF vectors are identical or not. As all vectors have been normalized, they will be identical if the directions of vectors are the same. The cosine of the angle between $\tilde{\mathbf{X}}\left(u,v\right)$ and $\tilde{\mathbf{X}}\left(\psi ,\omega \right)$ can be calculated by

$$\begin{array}{ll}\cos \left[\tilde{\mathbf{X}}\left(u,v\right),\hspace{0.17em}\tilde{\mathbf{X}}\left(\psi ,\omega \right)\right]& =\frac{{\tilde{\mathbf{X}}}^{\mathrm{T}}\left(u,v\right)\hspace{0.17em}\tilde{\mathbf{X}}\left(\psi ,\omega \right)}{{\Vert \tilde{\mathbf{X}}\left(u,v\right)\Vert }_2{\Vert \tilde{\mathbf{X}}\left(\psi ,\omega \right)\Vert }_2}\\ & =\langle \tilde{\mathbf{X}}\left(u,v\right),\hspace{0.17em}\tilde{\mathbf{X}}\left(\psi ,\omega \right)\rangle \end{array},$$

(9)

where $\langle \tilde{\mathbf{X}}\left(u,v\right),\hspace{0.17em}\tilde{\mathbf{X}}\left(\psi ,\omega \right)\rangle$ is the scalar product of $\tilde{\mathbf{X}}\left(u,v\right)$ and $\tilde{\mathbf{X}}\left(\psi ,\omega \right)$.

Therefore, Equation (8) will hold if

$$1-\langle \tilde{\mathbf{X}}\left(u,v\right),\hspace{0.17em}\tilde{\mathbf{X}}\left(\psi ,\omega \right)\rangle =0$$

(10)

Noise effect is not considered in the above derivation. In noisy environments, we cannot find SSPs that exactly satisfies Equation (10). Instead, we can get SSPs from the following criterion:

$$\left|1-\langle \tilde{\mathbf{X}}\left(u,v\right),\hspace{0.17em}\tilde{\mathbf{X}}\left(\psi ,\omega \right)\rangle \right|<{\delta }_1,$$

(11)

where $|·|$ is the absolute value of $·$, and ${\delta }_1$ is the SSPs threshold close to zero. Therefore, both $\tilde{\mathbf{X}}\left(u,v\right)$ and $\tilde{\mathbf{X}}\left(\psi ,\omega \right)$ are regarded as SSPs if they satisfy Equation (11).

As stated in [20], most of the signal energy will be concentrated in nearly 10% of the frequency bins. Therefore, in our study, the frequency bins are sorted in descending order according to their variances and only $N_f$ frequency bins with larger variance are selected for identifying SSPs. Moreover, we recommend that the data could be segmented when the sampling length is very large, and the results obtained in different segments can be combined by the similarity of the signal itself.

SSPs threshold ${\delta }_1$ has a large effect on the accuracy of SSPs and we now discuss how to choose the SSPs threshold. If ${\delta }_1$ is too small, the accuracy of SSPs will increase, however, the number of identified SSPs each time will decrease so that the efficiency of identifying enough SSPs will decline. Too small thresholds will even lead to insufficient SSPs. Otherwise, if ${\delta }_1$ is too large, the criterion becomes loose and too many outliers will be misjudged as SSPs, which will reduce the accuracy of mixing matrix estimation and source recovery. Since ${\delta }_1$ is related to the property of source signals, it is hard to give a unified range for all kinds of signals. A feasible approach that considers both efficiency and accuracy is to set a smaller threshold ${\delta }_1$ and a minimum number $N_{\mathrm{min-SSPs}}$ of identified SSPs. If the number $N_{\mathrm{SSPs}}$ of extracted SSPs is smaller than $N_{\mathrm{min-SSPs}}$, the threshold will be doubled. When source signals contain large noise or source signals are not very sparse in the TF domain, the threshold will gradually increase, which can reduce the effect of using an unsuitable threshold.

In general, TF vectors with negligible energy are greatly influenced by noise, which will easily lead to the misjudgment of SSPs. To obtain effective SSPs, these vectors should be removed before identifying SSPs if they satisfy

$${\Vert \mathbf{X}\left(u_0,v_0\right)\Vert }_2<{\delta }_2\overline{{\Vert \mathbf{X}\Vert }_2},$$

(12)

where ${\delta }_2$ is a threshold close to zero and $\overline{{\Vert \mathbf{X}\Vert }_2}$ represents the average of 2-norm of all TF vectors.

After identifying SSPs, the next stage is to estimate the mixing matrix. It can be found from Equations (6) and (7) that the identified SSPs are the set of normalized column vectors of the mixing matrix. Therefore, the mixing matrix can be estimated by clustering these TF vectors and the hierarchical clustering technique [21,22] is used here. It should be noted that this clustering technique may not be the best algorithm to cluster SSPs as other algorithms can also be used [23]. More details on adjusting the cluster number can be found in [17]. As studied in [17], the mixing matrix estimation error can be further reduced by removing the points which are away from the mean direction of the cluster. This strategy is also used in our study and the outlier detection rule is the same as [17]. By re-clustering SSPs after elimination of the outliers, each column of $\widehat{A}$ can be obtained via calculating the center of each cluster.

2.3. Source Recovery

Assumption 4.

At each TF point, the number of source signals is smaller than that of mixed signals.

Even though $\widehat{A}$ is known, the solution of the system in Equation (1) is not unique. Actually, source signals can be recovered by a series of least square problems [19] with Assumption 3 and Assumption 4, which minimizes the error function by selecting the optimal $N\times \left(N-1\right)$ submatrix of $\widehat{A}$. Let $\mathcal{A}$ be a set composed of all $N\times \left(N-1\right)$ submatrices of $\widehat{A}$, that is

$$\mathcal{A}=\left\{{\mathcal{A}}_{\hspace{0.17em}i}|{\mathcal{A}}_{\hspace{0.17em}i}=\left[{\widehat{a}}_{{\theta }_1},{\widehat{a}}_{{\theta }_2},\dots ,{\widehat{a}}_{{\theta }_{N-1}}\right]\right\}$$

(13)

Then for each TF point $\left(t,\hspace{0.17em}k\right)$, there exists ${\mathcal{A}}_{\hspace{0.17em}\ast }=\left[{\widehat{a}}_{{\varphi }_1},{\widehat{a}}_{{\varphi }_2},\dots ,{\widehat{a}}_{{\varphi }_{N-1}}\right]$ that satisfies

$$\mathbf{X}\left(t,\hspace{0.17em}k\right)={\mathcal{A}}_{\hspace{0.17em}\ast }{\mathcal{A}}_{\hspace{0.17em}\ast }^{†}\mathbf{X}\left(t,\hspace{0.17em}k\right),$$

(14)

where ${}^{†}$ is the pseudo-inverse of a matrix. Then, source signals can be estimated by

$${\widehat{S}}_j\left(t,\hspace{0.17em}k\right)=\{\begin{array}{ll}\begin{array}{l}{e}_i,\\ 0,\end{array}& \begin{array}{l}\mathrm{if}\hspace{0.17em}j={\varphi }_i\\ \mathrm{otherwise}\end{array},\end{array}$$

(15)

where $e={\left[e_1,\hspace{0.17em}{e}_2,\hspace{0.17em}\dots ,\hspace{0.17em}{e}_{N-1}\right]}^{\mathrm{T}}={\mathcal{A}}_{\hspace{0.17em}\ast }^{†}\mathbf{X}\left(t,\hspace{0.17em}k\right)$, and ${\mathcal{A}}_{\hspace{0.17em}\ast }$ can be obtained by

$${\mathcal{A}}_{\hspace{0.17em}\ast }=\arg \hspace{0.17em}\underset{{\mathcal{A}}_{\hspace{0.17em}i}\in \mathcal{A}}{\min }{\Vert \mathbf{X}\left(t,\hspace{0.17em}k\right)-{\mathcal{A}}_{\hspace{0.17em}i}{\mathcal{A}}_{\hspace{0.17em}i}^{†}\mathbf{X}\left(t,\hspace{0.17em}k\right)\Vert }_2$$

(16)

Finally, the time domain of the estimated source signals $\widehat{\mathbf{S}}\left(t\right)$ can be easily obtained by inverse STFT.

3. Proposed Source Contribution Estimation Method

Unlike the determined BSS, mixed signals in UBSS usually cannot be linearly represented by the estimated source signals, due to noise and separation errors, i.e., there exists a residual between mixed signals and the estimated source signals. Therefore, the i-th mixed signal $x_i$ can be expressed as

$$x_i=w_i^{\mathrm{T}}\widehat{\mathbf{S}}+z_i,$$

(17)

where $w_i={\left[w_{i1},\hspace{0.17em}{w}_{i2},\hspace{0.17em}\dots ,\hspace{0.17em}{w}_{iM}\right]}^{\mathrm{T}}$ represents the coefficients and $z_i$ represents the residual. It should be noted that $x_i={\left[x_i\left(1\right),x_i\left(2\right),\cdots ,x_i\left(\mathrm{T}\right)\right]}_{1\times T}$ represents the whole discrete sequence of the i-th mixed signal and $\widehat{\mathbf{S}}={\left[s_1^{\mathrm{T}},s_2^{\mathrm{T}},\dots ,s_M^{\mathrm{T}}\right]}^{\mathrm{T}}$ is also the whole discrete sequences of M estimated source signals, that is, $\widehat{\mathbf{S}}$ is a $M\times T$ matrix. $z_i={\left[z_i\left(1\right),z_i\left(2\right),\cdots ,z_i\left(T\right)\right]}_{1\times T}$ represents the i-th residual signal that has the same dimension as $x_i$. Then, a problem arises that how much $\widehat{\mathbf{S}}$ is contained in $x_i$. This problem can be addressed based on minimum redundant error energy, i.e.,

$$\begin{array}{ll}\underset{w_i}{\min }& {\Vert z_i\Vert }_2^2\end{array}$$

(18)

Therefore, the problem can be transformed into how to find the optimal coefficients in Equation (18).

Let $f\left(w_i\right)={\Vert z_i\Vert }_2^2$, it is easy to see that $f\left(w_i\right)$ is a continuous differential function. From Equation (17), we can obtain

$$\begin{array}{ll}f\left(w_i\right)& ={\Vert x_i-w_i^{\mathrm{T}}\widehat{\mathbf{S}}\Vert }_2^2\\ & =\left[x_i-w_i^{\mathrm{T}}\widehat{S}\right]{\left[x_i-w_i^{\mathrm{T}}\widehat{S}\right]}^{\mathrm{T}}\\ & =x_i{x}_i^{\mathrm{T}}-x_i{\widehat{\mathbf{S}}}^{\mathrm{T}}{w}_i-w_i^{\mathrm{T}}\widehat{\mathbf{S}}{x}_i^{\mathrm{T}}+w_i^{\mathrm{T}}\widehat{\mathbf{S}}{\widehat{\mathbf{S}}}^{\mathrm{T}}{w}_i\end{array}$$

(19)

The derivative of $f\left(w_i\right)$ with respect to $w_i$ can be calculated by

$${\nabla }_{w_i}\left[f\left(w_i\right)\right]=2\left[\widehat{\mathbf{S}}{\widehat{\mathbf{S}}}^{\mathrm{T}}{w}_i-\widehat{\mathbf{S}}{x}_i^{\mathrm{T}}\right]$$

(20)

Let ${\nabla }_{w_i}\left[f\left(w_i\right)\right]=0$, we will have

$$w_i={\left(\widehat{\mathbf{S}}{\widehat{\mathbf{S}}}^{\mathrm{T}}\right)}^{-1}\widehat{\mathbf{S}}{x}_i^{\mathrm{T}}$$

(21)

From Assumption 2, source signals are mutually independent, therefore, the estimated source signals are also approximately mutually independent, that is, the rank of $\widehat{\mathbf{S}}$ is M. Therefore, the rank of $\widehat{\mathbf{S}}{\widehat{\mathbf{S}}}^{\mathrm{T}}$ is also M, that is, $\widehat{\mathbf{S}}{\widehat{\mathbf{S}}}^{\mathrm{T}}$ is a matrix with full rank. Thus, $f\left(w_i\right)$ has only one stationary point. The Hessian matrix of $f\left(w_i\right)$ is

$${\nabla }_{w_i}^{2}f\left(w_i\right)=2\widehat{\mathbf{S}}{\widehat{\mathbf{S}}}^{\mathrm{T}}$$

(22)

The Hessian matrix of $f\left(w_i\right)$ is nonnegative definite. Therefore, from Equations (21) and (22), the minimum value of $f\left(w_i\right)$ can be obtained at $w_{i\ast }={\left(\widehat{\mathbf{S}}{\widehat{\mathbf{S}}}^{\mathrm{T}}\right)}^{-1}\widehat{\mathbf{S}}{x}_i^{\mathrm{T}}$.

Based on the above analysis, the optimal combination coefficients of $\widehat{\mathbf{S}}$ to all mixed signals can be obtained. Then, source contributions can be quantitatively estimated using $w_{\ast }$. The following equation can be obtained.

$$\begin{array}{llll}{x}_{i-j}=x_i-w_{i\ast ,\hspace{0.17em}j}{\widehat{s}}_j,& i=1,2,\dots ,N& \mathrm{and}& j=1,2,\dots ,M\end{array},$$

(23)

where $x_{i-j}$ represents the vector of $x_i$ that subtracts the contribution of ${\widehat{s}}_j$, and $w_{i\ast ,\hspace{0.17em}j}$ is the j-th element of $w_{i\ast }$. From Assumption 2, there is no cross energy among source signals, then the contribution $C_{ij}$ of the j-th estimated source signals ${\widehat{s}}_j$ to the i-th mixed signals $x_i$ are calculated according to

$$C_{ij}=1-\frac{{\Vert x_{i-j}\Vert }_2^2}{{\Vert x_i\Vert }_2^2}$$

(24)

Generally, due to the noise and estimation error, the sum of the contributions of all estimated source signals to a mixed signal is not equal to 1, which is different from complete BSS. From Equation (24), $C_{ij}<0$ if ${\Vert x_{i-j}\Vert }_2^2>{\Vert x_i\Vert }_2^2$, which implies that the j-th estimated source signals can decrease the overall vibration energy.

The flowchart of the proposed UBSS-based source contribution estimation method is shown in Figure 1.

4. Numerical Case Study

4.1. Performance of the Proposed UBSS Method

In this section, we mainly evaluate the separation performance of the proposed UBSS method with different sample sizes and different numbers of mixed signals. Some numerical case studies are conducted using five artificial source signals: $s_1\left(t\right)$ is a low frequency sinusoidal wave; $s_2\left(t\right)$ is a high frequency sinusoidal wave; $s_3\left(t\right)$ is a periodic wave with amplitude modulation; $s_4\left(t\right)$ is a shock attenuation signal wave; $s_5\left(t\right)$ is also a periodic wave with amplitude modulation. The generating functions of the source signals are listed as follows:

$$s(t)=\left[\begin{array}{l}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\\ s_4\left(t\right)\\ s_5\left(t\right)\end{array}\right]=\left[\begin{array}{l}\sin (2\pi \cdot 23t)\hfill \\ \sin (2\pi \cdot 281t)\hfill \\ \sin (2\pi \cdot 10t)\cos (2\pi \cdot 105t)\hfill \\ \begin{array}{l}5{\displaystyle \sum _{n=0}^{44}\left\{\begin{array}{l}\sin [2\pi \cdot 995t(t-0.023n)]\\ \exp [-628(t-0.023n)]u(t-0.023n)\end{array}\right\}}\\ \sin (2\pi \cdot 10t)\cos (2\pi \cdot 710t)\end{array}\hfill \end{array}\right]$$

(25)

Two, three and four mixed signals are generated by these five source signals. In each case, the averages of 50 Monte Carlo simulations are used to evaluate the performance of the proposed method, and in each simulation, Gaussian white noise with SNR = 10 dB is independently added to each source signal. The sampling frequency is 10 kHz. In the proposed method, the window length is 1024 and window overlap is 256, the number of selected frequency bins $N_f=80$, initial SSPs threshold δ₁ = 0.0001, minimum number of SSPs $N_{\mathrm{min-SSPs}}=300$ and energy threshold ${\delta }_2=0.1$.

To quantitatively verify the better performance of the proposed method, SNRs of $\widehat{A}$ and $\widehat{s}(t)$ are calculated by Equations (26) and (27), respectively.

$$E_A(i)=10\log \left({\Vert {\widehat{a}}_i\Vert }_2^2/{\Vert {\widehat{a}}_i-a_i\Vert }_2^2\right),$$

(26)

where $a_i$ and ${\widehat{a}}_i$ are the i-th column of $A$ and $\widehat{A}$, respectively.

$$E_s(i)=10\log \left(\underset{\varsigma }{\min }\hspace{0.17em}{\Vert {\widehat{s}}_i(t)\Vert }_2^2/{\Vert {\widehat{s}}_i(t)-\varsigma s_i(t)\Vert }_2^2\right),$$

(27)

where ς is a scalar that reflects the scalar indeterminacies.

The average SNRs of the estimated mixing matrix and the estimated source signals are shown in Figure 2a,b, respectively. From Figure 2a, the average SNRs of the estimated mixing matrix will increase with the increase in sample sizes. However, from Figure 2b, the average SNRs of the estimated source signals remained nearly unchanged with the increase in sample size. This is because the average SNRs of the estimated mixing matrix have been more than 40 dB when the sample size is 10,000, which means the mixing matrix is nearly the same as true mixing matrix. It can also be seen from Figure 2 that the separation performance also improves with the increase in the number of mixed signals. Though the average SNRs of the estimated mixing matrix with two mixtures are nearly the same as that with three mixtures, the average SNRs of the estimated source signals differ a lot in these two cases. That is because the number of source signals must be smaller than that of mixed signals at each TF point according to Assumption 4, which means that at most one source exists at each TF point in the case of two mixtures. This restriction is too strict, leading to worse separation performance in source signals.

4.2. Performance of the Proposed Source Contribution Estimation Method

In order to validate the effectiveness of the proposed source contribution quantitative estimation method, the following simulations are conducted. Source signals are the first four signals in Equation (25). The mixing matrix are

$$A=\left[\begin{array}{rrrr}\hfill 0.5736& \hfill 0.6293& \hfill 0.6561& \hfill 0.7071\\ \hfill 0.5587& \hfill -0.5953& \hfill -0.6536& \hfill 0.5930\\ \hfill 0.5991& \hfill 0.4995& \hfill -0.3774& \hfill -0.3851\end{array}\right]$$

(28)

The sampling frequency and sampling length is 10 kHz and 1 s, respectively. One hundred Monte Carlo simulations are conducted to evaluate the performance of the proposed method. In each simulation, Gaussian white noise is independently added into each source signal and each mixed signal with SNR = 10 dB and SNR = 15 dB, respectively.

The performance of the proposed UBSS method is compared with Reju’s method [17] and Zhen’s method [19]. Since Reju’s method is designed only for mixing matrix estimation, it cannot recover source signals. Therefore, the mixing matrix estimated by Reju’s method is then inputted into Zhen’s method to estimate source signals. The parameters in different methods are as follows. In all methods, the Hanning window is used in STFT, and the window length is 1024 and window overlap is 256. In Zhen’s method, regularization parameter $\lambda =0.001$ and energy threshold ${\delta }_2=0.1$. In Reju’s method, the parameter $\Delta \theta$ is set as ${1.5}^{\circ }$ and the number of selected frequency bins $N_f=80$. The parameters in the proposed method are the same as those in Section 4.1.

One example of the separation results is as follows. Waveforms and Fourier spectrums of source signals are displayed in Figure 3, while the major frequencies of the source signals can be easily obtained from Figure 3b. The major frequencies of $s_1(t)$, $s_2(t)$, and $s_4(t)$ are 23 Hz, 281 Hz, and 43 Hz, respectively, while the major frequencies of $s_3(t)$ are 95 Hz and 115 Hz. Waveforms and Fourier spectrums of mixed signals are shown in Figure 4. From Figure 4a, mixed signals are the superposition of source signals, therefore, we cannot directly obtain the waveforms of source signals. From Figure 4b, the major frequencies of source signals can be found in each Fourier spectrum of mixed signals, and the frequencies of $s_4(t)$ is overwhelmed by those of other source signals. Therefore, signal processing methods are needed to estimate all source signals.

The estimated mixing matrix of the proposed method is

$$\widehat{A}=\left[\begin{array}{rrrr}\hfill 0.5720& \hfill 0.6331& \hfill 0.6586& \hfill 0.6975\\ \hfill 0.5558& \hfill -0.5950& \hfill -0.6624& \hfill 0.6032\\ \hfill 0.6029& \hfill 0.4944& \hfill -0.3568& \hfill -0.3847\end{array}\right]$$

(29)

The absolute differences between $A$ and $\widehat{A}$ are calculated in Equation (30), illustrating that the mixing matrix has been well estimated because each of the absolute differences is very small.

$$\left|A-\widehat{A}\right|=\left[\begin{array}{rrrr}\hfill 0.0016& \hfill 0.0038& \hfill 0.0025& \hfill 0.0096\\ \hfill 0.0028& \hfill 0.0003& \hfill 0.0089& \hfill 0.0102\\ \hfill 0.0039& \hfill 0.0051& \hfill 0.0206& \hfill 0.0005\end{array}\right]$$

(30)

Source signals estimated by the proposed method, Zhen’s method and Reju’s method are displayed in Figure 5, Figure 6 and Figure 7, respectively. The order of $\widehat{s}(t)$ has been adjusted according to $s(t)$. Comparing Figure 5a with Figure 3a, we could find that the waveforms of $\widehat{s}(t)$ are quite similar to those of $s(t)$. From the Fourier spectrums of $\widehat{s}(t)$, the major frequencies of $s(t)$ have been well recovered, which can validate the effectiveness of the proposed UBSS method. As revealed by Figure 6a, it seems that the waveforms of $s(t)$ are also well recovered by Zhen’s method. However, as shown in Figure 6b, there is interference frequency 23 Hz in the Fourier spectrums of ${\widehat{s}}_3(t)$, and interference frequencies 95 Hz, 115 Hz and 281 Hz in the Fourier spectrums of ${\widehat{s}}_4(t)$, which indicates that source signal ${\widehat{s}}_3(t)$ and ${\widehat{s}}_4(t)$ were not well estimated. It could be seen from Figure 7 that $s_4(t)$ is not estimated by Reju’s method.

Average SNRs of 100 Monte Carlo simulations of $\widehat{A}$ estimated by different methods are listed in

Table 1. Comparison of average SNRs of estimated mixing matrix.

Methods	SNR (dB)				Average SNR of All Columns
	$a_1$	$a_2$	$a_3$	$a_4$
Zhen’s method	10.01	16.50	20.05	25.92	18.12
Reju’s method	39.63	38.47	25.35	26.15	32.40
The proposed method	43.65	41.82	37.44	39.69	40.65

, from which we can see that SNRs of $\widehat{A}$ estimated by the proposed method are larger than those estimated by Zhen’s method and Reju’s method. Average SNRs of all columns of the mixing matrix estimated by Zhen’s method, Reju’s method and the proposed method is 18.12 dB, 32.41 dB and 40.65 dB, respectively, which implies that the proposed method could estimate the mixing matrix more accurately.

Table 2. Comparison of average SNRs of estimated source signals.

Methods	SNR (dB)				Average SNR of All Sources
	${\widehat{\mathit{s}}}_1(\mathit{t})$	${\widehat{\mathit{s}}}_2(\mathit{t})$	${\widehat{\mathit{s}}}_3(\mathit{t})$	${\widehat{\mathit{s}}}_4(\mathit{t})$
Zhen’s method	9.61	9.65	7.39	6.99	8.41
Reju’s method	10.18	9.94	8.00	8.57	9.17
The proposed method	12.56	12.40	9.78	11.93	11.66

shows the average SNRs of 100 Monte Carlo simulations of $\widehat{s}\left(t\right)$ estimated by different methods. As can be seen in Table 2, all SNRs of $\widehat{s}\left(t\right)$ estimated by the proposed method are also larger than those estimated by Zhen’s method and Reju’s method. Average SNRs of all sources of Zhen’s method, Reju’s method and the proposed method are 8.41 dB, 9.17 dB and 11.66 dB, that is, the average SNR increments of all sources estimated by the proposed method are 38.72% and 27.18% when compared with Zhen’s method and Reju’s method, respectively. The above results tend to validate that the proposed UBSS method performs more effectively than Zhen’s method and Reju’s method.

The running time is used to evaluate the efficiency of the methods. CPU of the computer is Inter Core i5-4590 of 3.30 GHz and RAM is 1333 MHz DDR3 of 16 GB. Average time costs of the proposed method, Zhen’s method and Reju’s method are 1.79 s, 14.22 s and 0.17 s, respectively. The main difference between these three methods is the process of SSPs identification, which is the main cause for a significant difference in time cost. Reju’s method can identify SSPs according to single SSP, and only TF vectors in some frequency bins with a larger variance are selected for SSPs identification, therefore, time cost of Reju’s method is the least. SSPs must be identified between two TF vectors in Zhen’s method and the proposed method, which means more time consumption. However, SSPs are also identified in some frequency bins with a larger variance in the proposed method and they can be directly identified by searching the identical normalized TF vectors, instead of finding the sparsest coefficients. Therefore, the time cost of the proposed method is shorter than that of Zhen’s method.

Table 3. Average contribution comparison of estimated source signals.

Mixed Signals	Methods	Contributions (%)
		${\mathit{s}}_1(\mathit{t})$	${\mathit{s}}_2(\mathit{t})$	${\mathit{s}}_3(\mathit{t})$	${\mathit{s}}_4(\mathit{t})$
$x_1\left(t\right)$	Zhen’s method	0.2015	0.2790	0.1227	0.1049
	Reju’s method	0.2180	0.2790	0.1544	0.2379
	The proposed method	0.2368	0.2857	0.1734	0.3090
	Real contributions	0.2314	0.2780	0.1630	0.3018
$x_2\left(t\right)$	Zhen’s method	0.2152	0.2799	0.1534	0.1083
	Reju’s method	0.2591	0.2877	0.1665	0.2090
	The proposed method	0.2615	0.2963	0.1998	0.2517
	Real contributions	0.2529	0.2872	0.1862	0.2447
$x_3\left(t\right)$	Zhen’s method	0.3791	0.2791	0.1400	0.2591
	Reju’s method	0.4275	0.2810	0.0854	0.1239
	The proposed method	0.4381	0.3079	0.1009	0.1594
	Real contributions	0.4229	0.2939	0.0906	0.1501

shows the average results of source contributions quantitative estimation using different methods, including also the real contributions. It can be clearly seen that source contributions of the proposed method are closer to the real source contributions than those of Zhen’s method and Reju’s method. The average absolute errors of source contributions are also calculated and listed in

Table 4. Comparison of average contribution errors of estimated source signals.

Mixed Signals	Methods	Contribution Errors (%)
		${\mathit{s}}_1(\mathit{t})$	${\mathit{s}}_2(\mathit{t})$	${\mathit{s}}_3(\mathit{t})$	${\mathit{s}}_4(\mathit{t})$
$x_1\left(t\right)$	Zhen’s method	6.52	2.22	6.57	21.45
	Reju’s method	1.70	1.28	3.89	6.72
	The proposed method	0.84	0.78	1.13	1.37
$x_2\left(t\right)$	Zhen’s method	6.00	2.71	6.10	16.32
	Reju’s method	1.26	0.92	4.47	4.09
	The proposed method	1.05	1.00	1.37	1.15
$x_3\left(t\right)$	Zhen’s method	6.14	3.26	5.84	13.51
	Reju’s method	1.11	2.25	2.03	3.08
	The proposed method	1.80	1.78	1.03	1.04

. As revealed by the data in Table 4, most of the average contribution errors of the proposed method are the smallest among these three methods, implying that the proposed method has higher accuracy in source contribution. All contribution errors of the proposed method are less than 1.80%, however, three contribution errors are larger than 10% in Zhen’s method and three contribution errors are larger than 4% in Reju’s method. Actually, the accurate source estimation is the premise for correct contribution estimation. Therefore, it can be concluded that the proposed method performs more effectively in recovering source signals and quantitatively estimating source contributions.

5. Experimental Study with Cylindrical Structure

Some practical mechanical systems or their partial sections have the shape of cylindrical shells, such as the underwater vehicles. Generally, the sound radiation of underwater vehicles strongly interferes with their performance and safety. Therefore, it is quite important for underwater vehicles to reduce their radiation noise to accomplish tasks successfully. Before that, it is necessary to estimate sound sources in advance. When the number of sensors is smaller than that of sources, UBSS is an excellent method to estimate sources in these cases. Therefore, a test bed with a cylindrical shell structure is used to examine the effectiveness of the proposed method.

In the experiments, an adjustable speed motor is used as a vibration source and an eccentric mass disc is driven by the motor to simulate the unbalanced vibration. Two loudspeakers are also used to simulate two radiated noise sources and two arbitrary waveform generators are used to produce two different source signals which are the inputs of these two loudspeakers, respectively. Besides, mixed signals are collected by four sound pressure sensors and are recorded by GEN2i high-speed data recorder. Schematic diagram and photos of the test site are displayed in Figure 8 and Figure 9, respectively.

The motor is running at 1740 r/m. Inputs of two loudspeakers, denoted by loudspeaker 1 and loudspeaker 2, respectively, are sine waves of 713 Hz and 917 Hz, respectively. The sampling length and the sampling frequency are 10 s and 5000 Hz, respectively. The second and the fourth mixed signals are selected to estimate three source signals and only a section of data from 4 s to 6 s is used. Waveforms and Fourier spectrums of mixed signals are displayed in Figure 10. From Figure 10a, mixed signals are the superposition of source signals, therefore, we cannot directly obtain waveforms of source signals from mixed signals. From Figure 10b, the major frequencies of source signals can be found in each Fourier spectrum of mixed signals. Therefore, mixed signals need to be further processed to obtain source signals.

Source signals (displayed only from 4.5 s to 5 s) estimated by the proposed method, Zhen’s method and Reju’s method are illustrated in Figure 11, Figure 12 and Figure 13, respectively. As revealed in Figure 11b, the major frequencies of source signals estimated by the proposed method are 29 Hz, 937 Hz and 713 Hz, respectively, which are consistent with the frequencies set in the experiment. However, from Figure 12, both the major frequency 29 Hz of the motor and the major frequency 713 Hz of the loudspeaker 1 is mis-estimated in the same signal, as shown in the Fourier spectrum of the first separated signal. Therefore, 29 Hz and 713 Hz will be mistaken for coming from the same source using Zhen’s method. The major frequencies of the first signal estimated by Reju’s method are also 29 Hz and 713 Hz, as shown in Figure 13. The results tend to illustrate that source signals have been well estimated by the proposed method.

Reju’s method can identify SSP based on the character of single SSP, and the performance of this method will be degraded in noisy cases. In Zhen’s method, to increase the computational efficiency, SSPs are identified only in some TF vectors that are randomly selected from TF vectors of mixed signals. If no or very few SSPs corresponding to a source are included in selected TF vectors, this source will be estimated with large error. And this may be the main reason why the performance of Zhen’s method is not so good as that of the proposed method.

After obtaining the estimated source signals, their contributions to the mixed signals can be calculated and are presented in

Table 5. Contribution comparison of estimated source signals.

Mixed Signals	Methods	Contributions (%)
		${\mathit{s}}_1(\mathit{t})$	${\mathit{s}}_2(\mathit{t})$	${\mathit{s}}_3(\mathit{t})$
$x_1\left(t\right)$	Zhen’s method	21.82	46.50	39.00
	Reju’s method	28.29	43.30	12.98
	The proposed method	3.10	47.43	50.27
	Real contributions	7.31	47.15	47.64
$x_2\left(t\right)$	Zhen’s method	59.90	17.18	24.36
	Reju’s method	61.91	16.56	2.19
	The proposed method	38.97	16.96	40.20
	Real contributions	45.41	17.09	36.54

. Real source contributions also need to be obtained by the experiment [1]. When one source is stopped, the decreased amount in vibration energy of mixed signals is observed by the sensors. The decreased amount is regarded as the source contribution of the stopped source. The real source contributions are also given in Table 5. We can find the source contributions of the proposed method are closer to the real source contributions than those of Zhen’s method and Reju’s method.

The absolute errors of source contributions are calculated and listed in

Table 6. Comparison of contribution errors of estimated source signals.

Mixed Signals	Methods	Contribution Errors (%)
		${\mathit{s}}_1(\mathit{t})$	${\mathit{s}}_2(\mathit{t})$	${\mathit{s}}_3(\mathit{t})$
$x_1\left(t\right)$	Zhen’s method	14.51	0.65	8.64
	Reju’s method	20.98	3.85	34.66
	The proposed method	4.21	0.28	2.63
$x_2\left(t\right)$	Zhen’s method	14.49	0.09	12.18
	Reju’s method	16.50	0.53	34.35
	The proposed method	6.44	0.13	3.66

. All contribution errors of the proposed method are smaller than those of Zhen’s method and Reju’s method, implying that the proposed method has higher accuracy in source contribution. The largest contribution error of the proposed method is only 6.44%, however, three contribution errors are larger than 12% in Zhen’s method and four contribution errors are larger than 15% in Reju’s method. Actually, accurate estimation of source signals is the precondition for accurate estimation of computational complexity. To some extent, estimation accuracy of contribution increases with the increase in source signal estimation accuracy. As shown in Figure 12 and Figure 13, the first separated signals of Zhen’s method and Reju’s method contains major frequencies of two sources (motor and loudspeaker one). Therefore, the contributions of their first separated signal contain contributions of two real sources, which will lead to an increase in their source contributions. Besides, since part of the contributions of loudspeaker 1 is mis-assigned to their first separated signals, the contribution of their third separated signal will be smaller than the real contributions. Based on a more accurate estimation of source signals, the contribution errors of the proposed method are smaller than those of contrast methods.

Running time is also used to evaluate the efficiency of the methods. Time costs of the proposed method, Zhen’s method and Reju’s method are 1.86 s, 6.42 s and 0.32 s, respectively. Reju’s method could identify SSPs by the property of single SSP, which can have higher efficiency than Zhen’s method and the proposed method. The proposed method can identify SSPs only at some optimal frequency bins and SSPs are identified by directly searching the identical TF vectors in the selected frequency bins, which could the reason why the efficiency of the proposed method is higher than that of Zhen’s method.

It should be noted that the proposed method is designed for off-line processing system because it needs some data to identify SSPs. However, for a real-time monitoring system, we can process the data in a piecewise way by the proposed method. From the experiment, the running time of the proposed method is only 1.86 s to process the data with a length of 2 s. Therefore, we can split the data into a fixed length segment and analyze them by the proposed method.

After source signals are well recovered and source contributions are calculated, the influences of sources on mixed signals can be determined. The vibration sources estimated by the proposed method can be used to machinery condition monitoring and fault diagnosis when source signals are difficult to be directly obtained. The main vibration sources can also be determined according to their contributions. Therefore, some measures can be taken to reduce the impact of the main vibration sources.

6. Conclusions

To identify the major vibration and noise sources of the mechanical systems, a novel source contribution quantitative estimation method is proposed for UBSS. The accuracy of the source contribution results relies largely on the accuracy of source recovery. Only by recovering source signals more accurately can we obtain higher accuracy of source contribution estimation. From the results of numerical case studies, the proposed method can not only estimate source signals from their mixtures in underdetermined cases, but also quantitatively estimate the source contributions with average deviations <2%. The results of experimental studies with a cylindrical structure also show the effectiveness of the proposed method in sources restoration and quantitative contribution estimation. The comparative results tend to validate that the proposed method performs more powerfully. Therefore, our method could serve as an effective and promising tool to investigate the major vibration and noise sources, which can benefit vibration and noise monitoring and control of mechanical systems.