Non-orthogonal joint block diagonalization based on the LU or QR factorizations for convolutive blind source separation

1. Introduction

Blind source separation (BSS) deals with the problem of finding both the unknown input sources and unknown mixing system from only observed output mixtures. BSS has recently become the focus of intensive research work due to its high potential in many applications such as antenna processing, speech processing and pattern recognition [1-3]. The recent successes of BSS might be also used in mechanical engineering [4-7]. In these reported applications for tackling the BSS, there were two kinds of BSS models i.e., instantaneous and convolutive mixture models. The instantaneous mixture models with a simple structure have been described in many papers and books [8-10]. However, when it came to deal with convolutive mixture signals, BSS might face a number of difficulties which seriously hindered its feasibility [11, 12]. There is currently an endeavor of research in separating convolutive mixture signals, yet no fully satisfying algorithms have been proposed so far.

Many approaches have been proposed to solve the convolutive BSS (CBSS) problem in recent years. One kind of prevailing and representative approaches is joint block diagonalization (JBD) which can produce potentially more elegant solution for CBSS in time domain [13]. In the article, we focus on the JBD problem, which has been firstly treated in [14] for a set of positive definite symmetric matrices. And then the same conditions also mentioned in [15]. To solve the JBD problem, Belouchrani have sketched several Jacobi strategies in [16-18]: the JBD problem was turn into a minimization problem which was processed by iterative methods; as a product of Givens rotations, each rotation made a block-diagonality criterion minimum around a fixed axis. Févotte and Theis [19, 20] have pointed out that the behavior of Jacobi approach was very much dependent on the initialization of the orthogonal basis and also on the choice of the successive rotations. Then they proposed some strategies to improve the efficiency of JBD. But there are also several critical constraints: the joint block-diagonalizer is an orthogonal (unitary in the complex case) matrix; the spatial pre-whitening is likely to lead to a larger error and, moreover, this error is unable to correct in the subsequent analysis. In [21] a gradient-based JBD algorithm have been used to achieve the same task but for non-unitary joint block-diagonalizers. This approach suffers from slow convergence rate since the iteration criteria possesses a fixed step size. To improve this shortcoming, some gradient-based JBD algorithms with optimal step size have been provided and studied in [22-24]. However, these algorithms are apt to converge to a local minimum and have low computational efficiency. To eliminate the degenerate solutions of the nonunitary JBD algorithm, Zhang [25] optimized a penalty term based weighted least-squares criterion. In [26], a novel tri-quadratic cost function was introduced, furthermore, an efficient algebraic method based on triple iterations has been used to search the minimum point of the cost function. Unfortunately, this method exists redundancy value and arises error when the mixture matrix is inversed. Some new Jacobi-like algorithms [27, 28] for non-orthogonal joint diagonalization have been proposed, but unfortunately cannot be used to solve the problem of block diagonalization.

Our purpose, here, is not only to tackle the problem of the approximate JBD by discard the orthogonal constraint on the joint block-diagonalizer i.e., impose as few assumptions as possible on the matrix set, but also to propose the JBD algorithms characterized by the merits of simplicity, effectiveness, and computational efficiency. Subsequently, we suggest two novels non-orthogonal JBD algorithms as well as the Jacobi-type schemes. The new methods are an extension of joint diagonalization (JD) algorithms [29] based on LU and QR decompositions mentioned in [30, 31] to the block diagonal case.

This article is organized as follows: the JBD problem is stated in Section 2.1. The two proposed algorithms are derived in section 2.2, whose convergence is proved in Section 2.3. In Section 3, we show how to apply the JBD algorithms to tackle the CBSS problem. Section 4 gives the results of numerical simulation by comparing the proposed algorithms with the state-of-the-art gradient-based algorithms introduced in [23]. In Section 5, the novel JBD algorithms are proved to be effective in separating vibration source of convolutive type, which outperforms JBD_OG and JRJBD algorithms [19].

2. Problem formulation and non-orthogonal JBD algorithms

3. Application to CBSS

The problems of convolutive BSS (CBSS) occur in various applications. One typical application is in blind separation of vibration signals, which is fully studied in this paper for detecting the solution of the CBSS problems. The CBSS consists of estimating a set of unobserved source signals from their convolutive mixtures without requiring a priori knowledge of the sources and corresponding mixing system. Then the CBSS can be identified by means of JBD of a set of covariance matrices. We consider the following discrete-time MIMO model [25]:

where $t$ is the discrete time index, $l^{\text{'}}=1,\dots ,L^{\text{'}}\text{,}$$L^{\text{'}}$denotes FIR filter’s length. $\mathbf{S}\left(t\right)={\left[s_1\left(t\right),\cdot \cdot \cdot ,s_m\left(t\right)\right]}^T$ denotes source signal vector with the source numbers are $m$, and $\mathbf{X}\left(t\right)={\left[x_1\left(t\right),\cdot \cdot \cdot ,x_n\left(t\right)\right]}^T$ is the mixing signal vector obtained by$n$observation signals. In the mixing linear time-invariant system, the matrix-type impulse response $\mathbf{H}\left(l^{\text{'}}\right)=\left[{\mathbf{h}}_1\left(l^{\text{'}}\right),\cdot \cdot \cdot ,h_m\left(l^{\text{'}}\right)\right]$ consists of channel impulse responses $h_{ij}\left(l^{\text{'}}\right)$$(i=1,\dots ,n\text{,}$$j=1,\dots ,m)$. Aiming to the received signal on the $i$th array element, we take the $P+1$ sliding window and constitute a column vector:

Then putting the array element processed by sliding window together and defining the observed signal vector:

Hence:

where $\stackrel{-}{\mathbf{S}}\left(t\right)=\left[{\mathbf{s}}_1\right(t),{\mathbf{s}}_2(t),\cdot \cdot \cdot ,{\mathbf{s}}_m(t){]}^T$, and:

${\mathbf{H}}_{ij}$ is block element of $\mathbf{H}$ which matrix dimension is $(P+1)\times (L^{\text{'}}+P)$, and:

The following assumptions concerning the above mixture model Eq. (14) have to be made to ensure that it is possible to apply the proposed algorithms to CBSS [25].

Assumption 1. The source signals are zero mean, individually correlated in time but mutually uncorrelated, i.e., for all time delay $\tau$, the cross correlation functions $r_{s_i{s}_j}(t,\tau )=0$, $\forall i\ne j$, and the auto-correlation functions $r_{s_i{s}_i}(t,\tau )\ne 0$, $i=1,\dots ,m$. Meanwhile, the source signals have a different auto-correlation function.

Assumption 2. The sensor noises $\mathbf{N}\left(t\right)$ are zero mean, independent identically distributed with the same power ${\sigma }_N^2$. The noises are assumed to be independent with the sources.

Assumption 3. The mixing matrix $\mathbf{H}$ is assumed as column full rank. This requires that the length of the sliding window $P$ satisfies $n(P+1)\ge m(P+L^{\text{'}})$.

Assumption 1 is the core assumption. As is shown in [25], this assumption enables us to separate the sources from their convolutive mixtures by diagonalizing the second-order statistic of the reformulated observed signals (this will be addressed below). Assumption 2 enables us to easily deal with the noise and Assumption 3 guarantees that the mixing system is invertible, therefore it is a necessary condition that the source signals can be completely separated.

Under these assumptions, the spatial covariance matrices of the observations satisfy:

where ${\tau }_k\in \left[\mathrm{0,1},2,\dots ,q-1\right](k=1,\dots ,K)$ is the successive time delays, ${\mathbf{R}}_N(t,t-{\tau }_k)$ is computed according to [26]. It can be deduced from the above assumptions that in Eq. (15) the matrices take the following forms, respectively:

where the block matrices in ${\mathbf{R}}_{\stackrel{-}{x}}\left({\tau }_k\right)$ and ${\mathbf{R}}_{\stackrel{-}{s}}\left({\tau }_k\right)$ have the following form:

where $r(-a,-b)=E\left(x_i\right(t-a\left)x_j^T\right(t-b\left)\right)$, ${\mathbf{R}}_{s_i{s}_i}\left({\tau }_k\right)$ have the similar form, which is the $(L^{\text{'}}+P)\times (L^{\text{'}}+P)$ matrix. According to the Eq. (15), a group of matrices $\mathbf{R}=\left\{{\mathbf{R}}_k,k=1,\dots ,K\right\}$ which can be block diagonalized, and satisfy ${\mathbf{R}}_k=\mathbf{A}{\mathbf{D}}_k{\mathbf{A}}^T$, ${\mathbf{D}}_k$ have diagonalization structure. Hence, The JBD method mentioned in section 2.2 can be used to solve CBSS problem. Once the joint block diagonalizer $\mathbf{B}$ is determined, the recovered signals are obtained up to permutation and a filter by:

It is worth mentioning that the indeterminacies of amplitude and permutation exist in JBD algorithms correspond to the well-known indeterminacies in CBSS. The correlation matrices ${\mathbf{R}}_{\stackrel{-}{x}}\left({\tau }_k\right)$ is actually replaced by their discrete time series estimate. To acquire a good estimate of the discrete correlation matrices, we may divide the observed sequences (the output of the reformulated model (15)) into the appropriate length of the sample.

4. Numerical simulations

Simulations are now provided to illustrate the behavior and the performance of the proposed JBD algorithms (LUJBD, QRJBD). We will also compare the proposed algorithms with the JBD_OG, JBD_ORG in the robustness and efficiency by generating random dates. To achieve these purposes, a set of real block-diagonal matrices ${\mathbf{D}}_k$ (for all $k=1,\dots ,K$) are devised from random entries with a Gaussian distribution of zero mean and unit variance. Then, random noise entries with a Gaussian distribution of zero mean and variance ${\sigma }_N^2$ will be added on the off-diagonal blocks of the previous matrices ${\mathbf{D}}_k$. A signal to noise ratio can be defined as $SNR=10\mathrm{l}\mathrm{o}\mathrm{g}(1/{\sigma }_N^2)$.

To measure the quality of the different algorithms, the following performance index is used [22]:

where ${\mathbf{G}}_{ij}$ ($i,j\in \left\{1,\cdot \cdot \cdot ,r\right\}$) denotes the $\left(i,j\right)$th block matrix of $\mathbf{G}=\tilde{\mathbf{B}}\mathbf{A}$. This index will be used in the CBSS, which can take into account the inherent indetermination of the BSS problem. It is clearly that the smaller the index performance $I_{conv}\left(\mathbf{G}\right)$, the better the separation quality. Regarding to the charts, $I_{conv}\left(\mathbf{G}\right)$ is often given in dB i.e., $I_{conv}\left(\mathbf{G}\right)dB=10\mathrm{l}\mathrm{o}\mathrm{g}\left(I_{conv}\right(\mathbf{G}\left)\right)$. In all the simulations, the true mixing matrix $\mathbf{A}$ or $\mathbf{H}$ has been randomly chosen with zero mean and unit variance.

In Fig. 1 and Fig. 2, we focus on the exactly determined case $M=$ 9, $L=$ 3, $K=$ 20, $SNR=$ 60 dB, the results have been averaged over 100 runs. Fig. 1 represents the percentage of successful runs, where a run is declared successful w.r.t. the following criteria satisfy: $I_{conv}<10^{-7}$, (st1), (st2). Fig. 2 represents the average running time per successful convergence. Comparing the LUJBD and QRJBD approaches with the state-of-the-art JBD_OG, and JBD_ORG., we can confirm that the approaches proposed in Section 2.2 improve the performance of JBD better than the gradient-based methods: the LUJBD and QRJBD methods converge to the global minimum more frequently, see Fig. 1, and faster see Fig. 2 than JBD_OG, JBD_ORG. Under the same terminated criteria, it can be observed that the LUJBD and QRJBD methods also outperform the JBD_OG, JBD_ORG methods, and the LUJBD method show the best performance. We can also conclude that the sensitivity of the different convergence termination criteria for different JBD methods is diverse and, moreover, the percentage of successful runs and average running time are also varying. In other words, we should choose appropriate terminate criteria which is able to obtain the accuracy of block diagonalization, goodness of success rate and convergence speed.

In Fig. 3 we focus on the exactly determined case $M=$ 9, $L=$ 3, $K=$ 20, $SNR=$ 60 dB, the results have been averaged over 100 runs. Because various approaches have different iteration time of each step, Here, we consider all methods converge to the same time, which is more reasonable than converge to certain iteration steps mentioned in [23, 25]. The evolution of the performance index versus the convergence time shows that the convergence performance of the LUJBD and QRJBD methods is better than the JBDOG and JBDORG method. The LUJBD and QRJBD algorithms cost less time when performance index $I_{conv}$ reaches a stable convergence level, and have smaller value of performance index $I_{conv}$ when all algorithms converge to same time. In other words, the BSS methods proposed in this paper possess less convergence time, higher precision of convergence, faster convergence speed.

Fig. 1The percentage of successful runs

Fig. 2The average running time per successful convergence

Fig. 3Average performance index Iconv versus the convergence time for LUJBD, QLJBD, JBDOG, JBDORG algorithms. SNR=60 and K= 20 matrices

In Fig. 4, we discuss the number of the matrices how to affect the performance. The results have been averaged over 100 runs. We devise the same stop criteria i.e., one of the terminate criteria (st1), (st2) and (st3) is satisfied. We set $M=$ 9, $L=$ 3, $SNR=$ 60 dB. The following observations can be made: the more matrices to be joint block-diagonalized, the better performance we can obtain. But the computational cost also increases when the number of matrix rises. Therefore, the choice of matrix number should combine the accuracy of JBD algorithms with complexity of JBD algorithms. From Fig. 4, the matrix number 20 is a better choice. The LUJBD algorithm with better convergence turned out to be slightly superior to other JBD algorithms.

Fig. 4Average Iconv versus the number of matrices for LUJBD, QLJBD, JBDOG, JBDORG algorithms. SNR= 60 and avarage 100 runs

Fig. 5Average Iconv versus SNR for LUJBD, QLJBD, JBDOG, JBDORG algorithms. K= 20 and avarage 100 runs

With the same stop criteria and other assumptions in Fig. 4, one can observe from Fig. 5 that when the SNR grows, the average performance of each algorithm becomes better except few fluctuation points. The noise sensitivity of LUJBD is slightly higher that the remaining three kinds of methods, however, for a given value of SNR, the average performance index of LUJBD and QLJBD is always better than that of two gradient-based methods.

Finally, in Fig. 6 and Fig. 7, we study separation performance versus matrix dimension $M$ and the block dimension $L$ for both algorithms ($K=$ 20, $SNR=$ 60 dB), the results have been averaged over 100 runs. One can observe that in the same block dimension $L$ case, the larger the matrix dimension $M$, the weaker the estimation accuracy of mixing matrix. And in the same matrix dimension $M$ case, the larger the block dimension $L$, the better the estimation accuracy of mixing matrix. Therefore, we can improve the performance index by increasing the number of the matrix $K$ and the block dimension $L$ when the dimension of target matrix increases.

Fig. 6Average Iconv versus the convergence time with different M and L for LUJBD algorithms. SNR=60 and K= 20 matrices

Fig. 7Average Iconv versus the convergence time with different M and L for QLJBD algorithms. SNR=60 and K= 20 matrices

5. Applying CBSS to the vibration source separation

The experimental model is a double-stiffened cylindrical shell depicted in Fig. 8, which is used to simulate the cabin model. Underwater vibration tests of the double-stiffened cylindrical shell were carried out in anechoic water pool with a length of 16 meters, a width of 8 meters and a height of 8 meters. In the double-stiffened cylindrical shell, three exciters were arranged in the front part (No. 1 excitation source), middle part (No. 2 excitation source), rear part (No. 3 excitation source), respectively, which were used to simulate vibration sources of the internal equipment. Twenty-nine vibration acceleration sensors were arranged in the inner shell, and four accelerometers containing abundant vibration information were arranged in the vicinity of each excitation point. The location of exciters and acceleration sensors were shown in Fig. 9. Only the vertical excitation and response were considered in this test, and the model was located underwater 3.75 m. During the test, three exciters were controlled on the shore, and each exciter was turned on separately or multiple exciters were operated simultaneously according to different test requirements. Three exciters launched a continuous sinusoidal signal with different excitation frequencies and same energy, the frequency was 5 kHz, 4 kHz, 3 kHz, corresponding to No. 1 exciter, No. 2 exciter and No. 3 exciter. The vibration data was collected when the exciter was in stable operation. The sampling frequency was 16384 Hz and the sampling time was 10 s.

Fig. 11(d), (e), (f) show the mixture signals obtained by three sensors (17, 20, 23) on the inner shell when all of the exciters act simultaneously. It is obvious that the mix-signals with mutual spectrum aliasing are not able to represent real vibration characteristics and be utilized directly. We can also demonstrate that a mixture of vibrations is most often of the convolutive type which is not prone to be tackled, and moreover, the independence among the vibration sources is often not satisfactory strictly. Therefore, it is difficult to separate mechanical vibration source using traditional source separation methods. However, we propose the Jacobi-type JBD algorithms based on the LU or QR factorizations in this article, which can overcome above shortcomings effectively.

To ensure the stability of the solution, the terminate criteria $\epsilon =abs(I_{off}^{t+t_1}-I_{off}^{t+t_1+1})<10^{-6}$$(t_1=0,\dots ,9)$ is selected. We select observed signals $n=$ 5 i.e., sensors 17, 20, 23, 26, 29，and source signals $m=$ 3. The model parameters including $L^{\text{'}}$, $P$, $K$ are selected to guarantee that the solution accuracy satisfy $I_{off}\le$ –35 dB. The filter length $L^{\text{'}}=$ 13, the sliding window $P=$ 17 and a set of covariance matrices $K=$ 30 with a time lag ${\tau }_k$ taking linearly spaced valued.

In Fig. 10, the evolution of the performance index $I_{off}$ versus the convergence time shows that the convergence performance of the LUJBD and QRJBD methods are superior to the JBDOG and JRJBD methods [20]. The LUJBD and QRJBD algorithms cost less time when $I_{off}<$ –20 dB, and have smaller performance index $I_{off}$ when the convergence time is same. In other words, the BSS methods proposed in this paper possess higher precision of convergence, faster convergence speed. The low computational accuracy and efficiency of the latter two algorithms are mainly due to following reasons: (1) the JBDOG algorithm generally suffers from slow convergence rate and is apt to converge to a local minimum; and the accuracy of blind source separation is often hindered by the inversion of ill-conditioned matrices. (2) the joint block-diagonalizer of JRJBD algorithm is an orthogonal matrix, the spatial pre-whitening which is likely to lead to a larger error need to be applied. Moreover, this error is unable to correct in the subsequent analysis.

Fig. 11(a), (b), (c) represents source signals acquired by the acceleration sensors near excitation point (Here, we choose acceleration sensors 1, 8, 14 which have a higher signal-to-noise ratio) when each of the exciter act respectively. Comparison between the recovered primary sources – see Fig. 11(g), (h), (i) for LUJBD model and (j), (k), (l) for QRJBD model – with the true one shows that the proposed separation algorithms only subjected to indeterminacies of the permutation and amplitude are effective.

Fig. 8The experimental cabin model

Fig. 9The location of exciters and acceleration sensors

Fig. 10Average performance index Iconv versus the convergence time for four algorithms. K= 30, L'= 13, P= 17

Fig. 11Separation results for known channel-Sources: a) s1, b) s2, c) s3. Mixtures: d) x1, e) x2, f) x3. Separation sources for LUJBD: g) s1', h) s2', i) s3'. Separation sources for QLJBD: j) s1'', k) s2'', l) s3''

a)

b)

c)

d)

e)

f)

g)

h)

i)

j)

k)

l)

6. Conclusions

In this article, to solve the convolutive BSS (CBSS) problem, we present a class of simple Jacobi-type JBD algorithms based on the LU or QR factorizations. Using Jacobi-type matrices we can replace high dimensional minimization problems with a sequence of simple one-dimensional problems. The two novel methods named LUJBD and QRJBD are no more necessarily orthogonal, positive definite or symmetric matrices. In addition, we propose a novel convergence criteria which can reflect relative change between off-block diagonal and block diagonal. And this criterion can be adapted to all of iterative methods for solving JBD problem, which can also give a more intuitive comparison between different methods. The computation of the optimal parameter requires solving a polynomial of degree 2 in the real domain in LUJBD and QRJBD methods, which is more effective than other JBD methods that need to solve polynomial of degree 4. Therefore, the two new algorithms will cost fewer times than other non-orthogonal JBD methods, and moreover, the convergence of these two algorithms is also guaranteed.

A series of comparisons of the proposed approaches with the state-of-the-art JBD approaches (JBD_OG and JBD_ORG) based on gradient algorithms are implemented by varieties of numerical simulations. The results show that the LUJBD and QRJBD methods converge to the global minimum more frequently, and faster than JBD_OG, JBD_ORG. Choosing appropriate terminate criteria is beneficial to obtain the accuracy of block diagonalization, goodness of success rate and convergence speed. It can be readily observed that the more target matrices selected to be joint block-diagonalized, the better performance we can obtain. But the computational cost also increases when the number of matrix rises. Therefore, the choice of matrix number should combine algorithm accuracy with complexity. We can also improve the performance index by increasing the block dimension and decreasing the matrix dimension. Then the two novel JBD algorithms and the JBD_OG, JRJBD methods for separating practical vibration sources are studied. We can conclude that the convergence performance and accuracy of the LUJBD and QRJBD methods are superior to the JBD_OG and JRJBD methods. Finally, Comparison the recovered primary sources with the true one demonstrates the validity of the proposed algorithms for separating the vibration signals of convolutive mixtures.