An Adaptive Algorithm for Calculating Crosstalk Error for Blind Source Separation

Abstract

The crosstalk error is widely used to evaluate the performance of blind source separation. However, it needs to know the global separation matrix in advance, and it is not robust. In order to solve these problems, a new adaptive algorithm for calculating crosstalk error is presented, which calculates the crosstalk error by a cost function of least squares criterion, and the robustness of the crosstalk error is improved by introducing the position information of the maximum value in the global separation matrix. Finally, the method is compared with the conventional RLS algorithms in terms of performance, robustness and convergence rate. Furthermore, its validity is verified by simulation experiments and real world signals experiments.

Keywords

blind source separation crosstalk error global separation matrix robustness

1 Introduction

Blind signal separation (BSS) technology originated from the famous “cocktail party” problem (Choi and Cichocki, 1997) and has been a key research issue since then (Bridwell et al., 2018; Yatabe and Kitamura, 2018). It attempts to separate the unknown original signals only by using the observations of their mixture without having any prior information of the source signals and channels. It is important to measure the effectiveness of the blind separation algorithms, since it affects not only the effect of the algorithms, but also the iterative process of the algorithms. Evaluating the quality of the waveforms of the signals and the similarity of the waveforms is the most straightforward way to judge the effectiveness of algorithm. However, this method can only be regarded as a qualitative analysis, which cannot quantitatively reflect the stability and convergence rate of the algorithm. Therefore, researchers have proposed some performance indexes, such as crosstalk error (Cichocki and Amari, 2002; Parra and Sajda, 2003; El-Sankary et al., 2015), correlation coefficient (Li et al., 2012) and signal-to-noise ratio (Laheld and Cardoso, 1993; Li and Sejnowski, 1995; Grellier and Comon, 1998; Nandakumar and Bijoy, 2014; Aroudi et al., 2016; Mirzal, 2017) etc. The consistence between the order of the source signals and the separated signals is needed for calculating most of the performance indexes, but the order of the separated signals is uncertain. The crosstalk error does not need to consider the correspondence, and it can reflect the separation effect directly.

The crosstalk error represents the approximate degree between the global separation matrix C and a diagonal matrix. When C is closer to a diagonal matrix or its permutation matrix, the smaller the crosstalk error, the more similar the separated signals are to the source signals. Theoretically, the crosstalk error should be a good performance index in the field of BSS (Macchi, 1993; Moreau and Macchi, 1996; Amari, 1997; Yang and Amari, 1997; Moreau and Macchi, 1998; Macchi and Moreau, 1999; Yang, 1999; Basak and Amari, 1999). However, it is very sensitive to numerical errors and can not be directly used to evaluate separation effects in many practical situations. On one hand, in order to calculate the crosstalk error, it needs to know C in advance, which means that it cannot directly be applied to real-world experiments. On the other hand, if the maximum value of each row of C is in the same column, the crosstalk error usually does not correctly evaluate the separation effect. In this case, small crosstalk errors do not indicate that the separation effect is good. Since each component of the separated signal is an approximate estimation of the same source signal component, the actual separation effect is poor. The above two shortcomings affect its practical application. In order to overcome the above problems, Li Zong have proposed a method to improve crosstalk error evaluation criterion by introducing correlation coefficients to measure the degree of correlation between the separated signal components and adding them to the calculation formula of crosstalk error (Li et al., 2012). Although this method improves the accuracy of the evaluation when the maximum values of two rows of C are in the same column, C is still unknown when calculating crosstalk error in the practical application. Zhang used the Least Mean Square (LMS) algorithm to estimate C in real time (Zhang et al., 2011). However, the convergence rate of the method is limited by the convergence step size of the algorithm, and once the step size is inappropriate, the algorithm will fail. At the same time, when the value changes subtly in C (Mansour et al., 2002), the crosstalk error may vary greatly. In other words, the crosstalk error is not robust.

In response to the above problems, a fully adaptive crosstalk error estimation model based on the recursive least square (RLS) algorithm is presented in this paper. This new approach consists of sifting the impact of step size on the estimation of crosstalk error. Furthermore, the robustness of the crosstalk error is improved by introducing the position information of the maximum value of each row in the global separation matrix C in the process of calculating the crosstalk error.

The paper is organized as follows. In Section 2, the shortcomings of the original performance index are analysed. The main contributions of the paper can be found in Section 3, where the fully adaptive crosstalk error estimation model is presented firstly, and then the position information of the maximum value of each row is introduced to improve the robustness of the crosstalk error. Finally, the whole process of the algorithm is explained. The validity of the proposed algorithm is verified by experimental analysis in Section 4. Finally, Section 5 contains some conclusions.

2 Problem Formulation

Suppose that an antenna array is made up of M units, therefore the instantaneous linear mixing model can be expressed as

\begin{matrix} x_{t} = A S_{t} + N_{t}, \end{matrix}

\[ {\mathbf{x}_{t}}=\mathbf{A}{\mathbf{S}_{t}}+{\mathbf{N}_{t}},\]

(1)

where

x_{t} \in R^{m}

${\mathbf{x}_{t}}\in {R^{m}}$

stands for the observation vector and

S_{t} = [s_{1} (t), s_{2} (t), \dots, s_{n} {(t)]}^{T}

${\mathbf{S}_{t}}=[{s_{1}}(t),{s_{2}}(t),\dots ,{s_{n}}{(t)]^{T}}$

is the source signal vector.

N_{t} = {[n_{1} (t), n_{2} (t), \dots, n_{m} (t)]}^{T}

${\mathbf{N}_{t}}={[{n_{1}}(t),{n_{2}}(t),\dots ,{n_{m}}(t)]^{T}}$

is the additive observation noise vector. The subscript in

x_{t}

${\mathbf{x}_{t}}$

S_{t}

${\mathbf{S}_{t}}$

and

N_{t}

${\mathbf{N}_{t}}$

stands for the vectors at the moment t.

A \in R^{m \times n}

$\mathbf{A}\in {R^{m\times n}}$

is named as mixed matrix. In general, the noise signal

n_{i} (t)

${n_{i}}(t)$

i = 1, 2, \dots, m

$i=1,2,\dots ,m$

satisfies Gaussian distribution, and the literature (Qian et al., 2006) pointed out that the separated signal after the blind source separation is only the superposition of source signal and Gaussian noise.

In order to simplify the problem of BSS and reduce the computational complexity, it is generally necessary to pre-whiten the signal, and the signal after whitening is

\begin{matrix} z_{t} = V x_{t}, \end{matrix}

\[ {\mathbf{z}_{t}}=\mathbf{V}{\mathbf{x}_{t}},\]

(2)

where V is a

n \times m

$n\times m$

whitening matrix. After whitening, the problem of BSS changes from seeking for the mixed matrix A to seeking an

n \times n

$n\times n$

orthogonal decomposition matrix W (Jiang, 2014). Supposing

y_{t}

${\mathbf{y}_{t}}$

is the estimated vector of

S_{t} = {[s_{1} (t), s_{2} (t), \dots, s_{n} (t)]}^{T}

${\mathbf{S}_{t}}={[{s_{1}}(t),{s_{2}}(t),\dots ,{s_{n}}(t)]^{T}}$

and

y_{t}

${\mathbf{y}_{t}}$

can be modelled as

\begin{matrix} y_{t} = W z_{t} = WVA S_{t} + WV N_{t} . \end{matrix}

\[ {\mathbf{y}_{t}}=\mathbf{W}{\mathbf{z}_{t}}=\mathbf{WVA}{\mathbf{S}_{t}}+\mathbf{WV}{\mathbf{N}_{t}}.\]

(3)

The global separation matrix C is defined as

\begin{matrix} C = WVA, \end{matrix}

\[ \mathbf{C}=\mathbf{WVA},\]

(4)

C is a generalized permutation matrix (Cichocki and Amari, 2002), which is equal to the product of a permutation matrix and a full rank diagonal matrix:

\begin{matrix} C = P Λ, \end{matrix}

\[ \mathbf{C}=\mathbf{P}\boldsymbol{\Lambda },\]

(5)

where P is a permutation matrix, which represents the uncertainty of the order of the blind separation results. Λ is a full rank diagonal matrix, which represents the uncertainty of the amplitude of the blind separation results. It can be seen that the waveforms of signal can be completely recovered except for the real amplitude and order of the source signals.

In practice, the blind source separation algorithms can only make C as close as possible to a generalized permutation matrix. Based on this, in order to evaluate the performance of the separation algorithms, the literature (Cichocki and Amari, 2002) gave a measurement method by using the difference between C and a generalized permutation matrix, which is named as crosstalk error and expressed as

\begin{matrix} PI (C) = \frac{1}{n (n - 1)} {\sum_{i = 1}^{n} (\sum_{j = 1}^{n} \frac{| c_{i j} |}{{max}_{k} | c_{i k} |} - 1) + \sum_{j = 1}^{n} (\sum_{i = 1}^{n} \frac{| c_{i j} |}{{max}_{k} | c_{k j} |} - 1)}, \end{matrix}

\[ \textit{PI}(\mathbf{C})=\frac{1}{n(n-1)}\Bigg\{{\sum \limits_{i=1}^{n}}\Bigg({\sum \limits_{j=1}^{n}}\frac{|{c_{ij}}|}{{\max _{k}}|{c_{ik}}|}-1\Bigg)+{\sum \limits_{j=1}^{n}}\Bigg({\sum \limits_{i=1}^{n}}\frac{|{c_{ij}}|}{{\max _{k}}|{c_{kj}}|}-1\Bigg)\Bigg\},\]

(6)

where

c_{i j}

${c_{ij}}$

stands for the element in row i and column j. The smaller the value of

PI

$\mathit{PI}$

, the better the separation effect of the algorithm. If and only if C is a generalized permutation matrix,

PI = 0

$\textit{PI}=0$

. According to formula (6), it can be seen that the key to compute the crosstalk error is to obtain the matrix C. Another shortcoming of the crosstalk error is that it is not robust. As an example, the following two global separation matrices are presented:

\begin{matrix} C_{1} = [\begin{array}{ccccc} 100000 & 10 & 1 & 2 & 4 \\ 1 & 0 & 0 & 0 & 0 \\ 10 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 8 & 0 & 0 & 0 & 0 \end{array}]; C_{2} = [\begin{array}{ccccc} 100000 & 10 & 1 & 2 & 4 \\ 1 & 0 & 0 & 0 & 0 \\ 10 & 0 & 0.01 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 8 & 0 & 0 & 0 & 0 \end{array}] . \end{matrix}

\[ {\mathbf{C}_{1}}=\left[\begin{array}{c@{\hskip4.0pt}c@{\hskip4.0pt}c@{\hskip4.0pt}c@{\hskip4.0pt}c}100000\hspace{1em}& 10\hspace{1em}& 1\hspace{1em}& 2\hspace{1em}& 4\\ {} 1\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\\ {} 10\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\\ {} 1\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\\ {} 8\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\end{array}\right];\hspace{1em}{\mathbf{C}_{2}}=\left[\begin{array}{c@{\hskip4.0pt}c@{\hskip4.0pt}c@{\hskip4.0pt}c@{\hskip4.0pt}c}100000\hspace{1em}& 10\hspace{1em}& 1\hspace{1em}& 2\hspace{1em}& 4\\ {} 1\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\\ {} 10\hspace{1em}& 0\hspace{1em}& 0.01\hspace{1em}& 0\hspace{1em}& 0\\ {} 1\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\\ {} 8\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\hspace{1em}& 0\end{array}\right].\]

(7)

In matrices

C_{1}

${\mathbf{C}_{1}}$

and

C_{2}

${\mathbf{C}_{2}}$

, the maximum value of every row is in the first column, so the isolated signal is only the same signal. The corresponding separation effect of

C_{1}

${\mathbf{C}_{1}}$

and

C_{2}

${\mathbf{C}_{2}}$

is very poor, and the value of

PI (C_{1})

$\mathit{PI}({\mathbf{C}_{1}})$

and

PI (C_{2})

$\mathit{PI}({\mathbf{C}_{2}})$

should be very big. However, the crosstalk errors of them are

PI (C_{1}) = 0.00002

$\mathit{PI}({\mathbf{C}_{1}})=0.00002$

and

PI (C_{2}) = 0.00057

$\mathit{PI}({\mathbf{C}_{2}})=0.00057$

, which are very small. At the same time, from matrices

PI (C_{1})

$\mathit{PI}({\mathbf{C}_{1}})$

and

PI (C_{2})

$\mathit{PI}({\mathbf{C}_{2}})$

, we can also see that the crosstalk error will change greatly even though there is a minor change in C.

Therefore, the maximum value of all the rows in the C should be in different columns. However, from the definition of formula (6), this restriction is not considered. Therefore, in this paper, we add the position information of the maximum value of each row in the C to the definition of the crosstalk error to improve the robustness of the crosstalk error.

3 A New Algorithm for Calculating the Crosstalk Error

For overcoming the shortcomings of the original crosstalk error, firstly a model for estimating the mixed matrix A is established and the corresponding analytic solution is given in this section. Then, the position information constraint is introduced into the definition of crosstalk error to improve its robustness. Finally, the whole process of the algorithm for calculating the crosstalk error is given.

3.1 Modelling and Model Solving

In this subsection, an estimation method based on recursive least-squares (RLS) is proposed to estimate the global separation matrix in real time. According to formula (4), if we get matrix A, then C can be calculated by using

C = WVA

$\mathbf{C}=\mathbf{WVA}$

. Since

y_{t} = WV x_{t}

${\mathbf{y}_{t}}=\mathbf{WV}{\mathbf{x}_{t}}$

, so we can get that

x_{t}

${\mathbf{x}_{t}}$

can also be regarded as a mixture of

y_{t}

${\mathbf{y}_{t}}$

. Suppose

\begin{matrix} x_{t} = \hat{A} y_{t}, \end{matrix}

\[ {\mathbf{x}_{t}}=\hat{\mathbf{A}}{\mathbf{y}_{t}},\]

(8)

where

\hat{A} \in R^{m \times n}

$\hat{\mathbf{A}}\in {R^{m\times n}}$

. According to formulas (3), (5) and (8), we can get

\begin{matrix} A = \hat{A} P Λ, \end{matrix}

\[ \mathbf{A}=\hat{\mathbf{A}}\mathbf{P}\boldsymbol{\Lambda },\]

(9)

\begin{matrix} C = WVA = WV (\hat{A} P Λ) = (WV \hat{A}) (P Λ) . \end{matrix}

\[ \mathbf{C}=\mathbf{WVA}=\mathbf{WV}(\hat{\mathbf{A}}\mathbf{P}\boldsymbol{\Lambda })=(\mathbf{WV}\hat{\mathbf{A}})(\mathbf{P}\boldsymbol{\Lambda }).\]

(10)

It can be seen that there are differences in amplitude and order between the mixed matrix A and its estimated value

\hat{A}

$\hat{\mathbf{A}}$

, which do not affect the result of the crosstalk error defined by formula (6). The following cost function for calculating

\hat{A}

$\hat{\mathbf{A}}$

is established by using the least squares criterion, which is

\begin{matrix} min_{{\hat{A}}_{t}} J ({\hat{A}}_{t}) = min_{{\hat{A}}_{t}} \sum_{i = 1}^{t} λ^{t - i} ‖ {\hat{A}}_{t} y_{t} - x_{t} ‖^{2}, \end{matrix}

\[ \underset{{\hat{\mathbf{A}}_{t}}}{\min }J({\hat{\mathbf{A}}_{t}})=\underset{{\hat{\mathbf{A}}_{t}}}{\min }{\sum \limits_{i=1}^{t}}{\lambda ^{t-i}}\| {\hat{\mathbf{A}}_{t}}{\mathbf{y}_{t}}-{\mathbf{x}_{t}}{\| ^{2}},\]

(11)

where λ is the forgetting factor subjecting to

0 < λ < 1

$0<\lambda <1$

and close to 1,

{\hat{A}}_{t}

${\hat{\mathbf{A}}_{t}}$

is the estimated value of A at moment t.

The gradient of

J ({\hat{A}}_{t})

$J({\hat{\mathbf{A}}_{t}})$

in (11) is

\begin{matrix} \frac{\partial J}{\partial {\hat{A}}_{t}} = - 2 \sum_{i = 1}^{t} λ^{t - i} {- x_{t} y_{t}^{T} + {\hat{A}}_{t} y_{t} y_{t}^{T}} . \end{matrix}

\[ \frac{\partial J}{\partial {\hat{\mathbf{A}}_{t}}}=-2{\sum \limits_{i=1}^{t}}{\lambda ^{t-i}}\big\{-{\mathbf{x}_{t}}{\mathbf{y}_{t}^{T}}+{\hat{\mathbf{A}}_{t}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}\big\}.\]

(12)

The zero point of formula (12) is the optimum solution of formula (11).

Let

\frac{\partial J}{\partial {\hat{A}}_{t}} = - 2 \sum_{i = 1}^{t} λ^{t - i} {- x_{t} y_{t}^{T} + {\hat{A}}_{t} y_{t} y_{t}^{T}} = 0

$\frac{\partial J}{\partial {\hat{\mathbf{A}}_{t}}}=-2{\textstyle\sum _{i=1}^{t}}{\lambda ^{t-i}}\{-{\mathbf{x}_{t}}{\textbf{y}_{t}^{T}}+{\hat{\mathbf{A}}_{t}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}\}=0$

, then we can get

\begin{matrix} {\hat{A}}_{t} = (\sum_{i = 1}^{t} λ^{t - i} x_{t} {y_{t}}^{T}) {(\sum_{i = 1}^{t} λ^{t - i} y_{t} y_{t}^{T})}^{- 1} . \end{matrix}

\[ {\hat{\mathbf{A}}_{t}}=\Bigg({\sum \limits_{i=1}^{t}}{\lambda ^{t-i}}{\mathbf{x}_{t}}{{\mathbf{y}_{t}}^{T}}\Bigg){\Bigg({\sum \limits_{i=1}^{t}}{\lambda ^{t-i}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}\Bigg)^{-1}}.\]

(13)

Since calculating the inverse matrix is complex, therefore an iterative method for calculating

{\hat{A}}_{t}

${\hat{\mathbf{A}}_{t}}$

is proposed. The update formula of

{\hat{A}}_{t}

${\hat{\mathbf{A}}_{t}}$

is discussed in the following.

Let

φ_{t} = \sum_{i = 1}^{t} λ^{t - i} y_{t} y_{t}^{T}

${\varphi _{t}}={\textstyle\sum _{i=1}^{t}}{\lambda ^{t-i}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}$

, then we can get the following formula:

\begin{matrix} φ_{t} = \sum_{i = 1}^{t - 1} λ^{t - 1 - i} y_{t - 1} y_{t - 1}^{T} + y_{t} y_{t}^{T} = λ φ_{t - 1} + y_{t} y_{t}^{T} . \end{matrix}

\[ {\varphi _{t}}={\sum \limits_{i=1}^{t-1}}{\lambda ^{t-1-i}}{\mathbf{y}_{t-1}}{\mathbf{y}_{t-1}^{T}}+{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}=\lambda {\varphi _{t-1}}+{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}.\]

(14)

Left multiplying

{φ_{t}}^{- 1}

${{\varphi _{t}}^{-\mathrm{1}}}$

and right multiplying

φ_{t - 1}^{- 1}

${\varphi _{t-1}^{-1}}$

on both sides of the equation (14), then

\begin{matrix} {φ_{t}}^{- 1} φ_{t} φ_{t - 1}^{- 1} = λ {φ_{t}}^{- 1} φ_{t - 1} φ_{t - 1}^{- 1} + {φ_{t}}^{- 1} y_{t} y_{t}^{T} φ_{t - 1}^{- 1}, \end{matrix}

\[ {{\varphi _{t}}^{-\mathrm{1}}}{\varphi _{t}}{\varphi _{t-1}^{-1}}=\lambda {{\varphi _{t}}^{-\mathrm{1}}}{\varphi _{t-1}}{\varphi _{t-1}^{-1}}+{{\varphi _{t}}^{-\mathrm{1}}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}{\varphi _{t-1}^{-1}},\]

(15)

\begin{matrix} φ_{t - 1}^{- 1} = λ {φ_{t}}^{- 1} + {φ_{t}}^{- 1} y_{t} y_{t}^{T} φ_{t - 1}^{- 1} . \end{matrix}

\[ {\varphi _{t-1}^{-1}}=\lambda {{\varphi _{t}}^{-\mathrm{1}}}+{{\varphi _{t}}^{-\mathrm{1}}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}{\varphi _{t-1}^{-1}}.\]

(16)

Right multiplying

y_{t}

${\mathbf{y}_{t}}$

on both sides of equation (16), then

\begin{matrix} φ_{t - 1}^{- 1} y_{t} = λ {φ_{t}}^{- 1} y_{t} + {φ_{t}}^{- 1} y_{t} y_{t}^{T} φ_{t - 1}^{- 1} y_{t} = {φ_{t}}^{- 1} y_{t} [λ + y_{t}^{T} φ_{t - 1}^{- 1} y_{t}], \end{matrix}

\[ {\varphi _{t-1}^{-1}}{\mathbf{y}_{t}}=\lambda {{\varphi _{t}}^{-\mathrm{1}}}{\mathbf{y}_{t}}+{{\varphi _{t}}^{-\mathrm{1}}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}{\varphi _{t-1}^{-1}}{\mathbf{y}_{t}}={{\varphi _{t}}^{-\mathrm{1}}}{\mathbf{y}_{t}}\big[\lambda +{\mathbf{y}_{t}^{T}}{\varphi _{t-1}^{-1}}{\mathbf{y}_{t}}\big],\]

(17)

\begin{matrix} {φ_{t}}^{- 1} y_{t} = \frac{φ_{t - 1}^{- 1} y_{t}}{[λ + y_{t}^{T} φ_{t - 1}^{- 1} y_{t}]} . \end{matrix}

\[ {{\varphi _{t}}^{-\mathrm{1}}}{\mathbf{y}_{t}}=\frac{{\varphi _{t-1}^{-1}}{\mathbf{y}_{t}}}{[\lambda +{\mathbf{y}_{t}^{T}}{\varphi _{t-1}^{-1}}{\mathbf{y}_{t}}]}.\]

(18)

Substituting equation (18) into (16), we can get

\begin{matrix} {φ_{t}}^{- 1} = \frac{1}{λ} [φ_{t - 1}^{- 1} - \frac{φ_{t - 1}^{- 1} y_{t} y_{t}^{T} φ_{t - 1}^{- 1}}{λ + y_{t}^{T} φ_{t - 1}^{- 1} y_{t}}] . \end{matrix}

\[ {{\varphi _{t}}^{-\mathrm{1}}}=\frac{1}{\lambda }\bigg[{\varphi _{t-1}^{-1}}-\frac{{\varphi _{t-1}^{-1}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}{\varphi _{t-1}^{-1}}}{\lambda +{\mathbf{y}_{t}^{T}}{\varphi _{t-1}^{-1}}{\mathbf{y}_{t}}}\bigg].\]

(19)

Let

P_{t} = φ_{t}^{- 1}

${\mathbf{P}_{\mathrm{t}}}={\varphi _{t}^{-1}}$

Q_{t} = P_{t - 1} / (λ + y_{t}^{T} P_{t - 1} y_{t})

${\mathbf{Q}_{t}}={\mathbf{P}_{t-1}}/(\lambda +{\mathbf{y}_{t}^{T}}{\mathbf{P}_{t-1}}{\mathbf{y}_{t}})$

q_{t} = \sum_{i = 1}^{t} λ^{t - i} x_{t} {y_{t}}^{T} = q_{t - 1} + x_{t} {y_{t}}^{T}

${\mathbf{q}_{t}}={\textstyle\sum \limits_{i=1}^{t}}{\lambda ^{t-i}}{\mathbf{x}_{t}}{{\mathbf{y}_{t}}^{T}}={\mathbf{q}_{t-1}}+{\mathbf{x}_{t}}{{\mathbf{y}_{t}}^{T}}$

, then it can be gotten that

\begin{matrix} P_{t} = \frac{1}{λ} [P_{t - 1} - P_{t - 1} y_{t} y_{t}^{T} Q_{t}] . \end{matrix}

\[ {\mathbf{P}_{t}}=\frac{1}{\lambda }\big[{\mathbf{P}_{t-1}}-{\mathbf{P}_{t-1}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}{\mathbf{Q}_{t}}\big].\]

(20)

Substitute

P_{t}

${\mathbf{P}_{t}}$

Q_{t}

${\mathbf{Q}_{t}}$

and

q_{t}

${\mathbf{q}_{t}}$

into (13), then

\begin{aligned} {\hat{A}}_{t} & = (\sum_{i = 1}^{t} λ^{t - i} x_{t} y_{t}^{T}) {(\sum_{i = 1}^{t} λ^{t - i} y_{t} y_{t}^{T})}^{- 1} \\ = \frac{1}{λ} (λ q_{t - 1} + x_{t} y_{t}^{T}) (P_{t - 1} - P_{t - 1} y_{t} y_{t}^{T} Q_{t}) \\ = q_{t - 1} P_{t - 1} - q_{t - 1} P_{t - 1} y_{t} y_{t}^{T} Q_{t} + x_{t} y_{t}^{T} P_{t} \\ = {\hat{A}}_{t - 1} - {\hat{A}}_{t - 1} y_{t} y_{t}^{T} Q_{t} + x_{t} y_{t}^{T} P_{t} . \end{aligned}

\[\begin{aligned}{}{\hat{\mathbf{A}}_{t}}& =\Bigg(\hspace{0.1667em}{\sum \limits_{i=1}^{t}}{\lambda ^{t-i}}{\mathbf{x}_{t}}{\mathbf{y}_{t}^{T}}\Bigg){\Bigg(\hspace{0.1667em}{\sum \limits_{i=1}^{t}}{\lambda ^{t-i}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}\Bigg)^{-1}}\\ {} & =\frac{1}{\lambda }\big(\lambda {\mathbf{q}_{t-1}}+{\mathbf{x}_{t}}{\mathbf{y}_{t}^{T}}\big)\big({\mathbf{P}_{t-1}}-{\mathbf{P}_{t-1}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}{\mathbf{Q}_{t}}\big)\\ {} & ={\mathbf{q}_{t-1}}{\mathbf{P}_{t-1}}-{\mathbf{q}_{t-1}}{\mathbf{P}_{t-1}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}{\mathbf{Q}_{t}}+{\mathbf{x}_{t}}{\mathbf{y}_{t}^{T}}{\mathbf{P}_{t}}\\ {} & ={\hat{\mathbf{A}}_{t-1}}-{\hat{\mathbf{A}}_{t-1}}{\mathbf{y}_{t}}{\mathbf{y}_{t}^{T}}{\mathbf{Q}_{t}}+{\mathbf{x}_{t}}{\mathbf{y}_{t}^{T}}{\mathbf{P}_{t}}.\end{aligned}\]

(21)

The global separation matrix C can be estimated by the following formula:

\begin{matrix} \hat{C} = WV {\hat{A}}_{t} . \end{matrix}

\[ \hat{\mathbf{C}}=\mathbf{WV}{\hat{\mathbf{A}}_{t}}.\]

(22)

An iterative method for calculating the crosstalk error can be obtained by substituting equation (22) into formula (6).

3.2 Location Introduction

It can be seen from the first section that some isolated signals may be the same signal, if the maximum value of the corresponding rows share the same column. Therefore, an improved crosstalk error definition is proposed here to eliminate the wrong decision caused by positions of optimal values in the same column, which is defined as

\begin{aligned} PI (C) & = \frac{1}{n (n - 1)} {\sum_{i = 1}^{n} (\sum_{j = 1}^{n} \frac{| c_{i j} |}{{max}_{k} | c_{i k} |} - 1) \\ + \sum_{j = 1}^{n} (\sum_{i = 1}^{n} \frac{| c_{i j} |}{{max}_{k} | c_{k j} |} \prod_{p = 1}^{j} (1 + r_{p}) - 1)}, \end{aligned}

\[\begin{aligned}{}\mathit{PI}(\mathbf{C})& =\frac{1}{n(n-1)}\Bigg\{{\sum \limits_{i=1}^{n}}\Bigg({\sum \limits_{j=1}^{n}}\frac{|{c_{ij}}|}{{\max _{k}}|{c_{ik}}|}-1\Bigg)\\ {} & \hspace{1em}+{\sum \limits_{j=1}^{n}}\Bigg({\sum \limits_{i=1}^{n}}\frac{|{c_{ij}}|}{{\max _{k}}|{c_{kj}}|}{\prod \limits_{p=1}^{j}}(1+{r_{p}})-1\Bigg)\Bigg\},\end{aligned}\]

(23)

where

r_{p} = 1

${r_{p}}=1$

, and

r_{p} = 0

${r_{p}}=0$

r_{p} > 0

${r_{p}}>0$

, when

p = 2, 3, \dots, n

$p=2,3,\dots ,n$

. Let

r_{p} > 0

${r_{p}}>0$

, if the column with the maximum value in row p is same as any other rows before row p, otherwise

r_{p} = 0

${r_{p}}=0$

, for

p = 2, 3, \dots, n

$p=2,3,\dots ,n$

r_{p}

${r_{p}}$

is a factor for punishing the proximity of C to a generalized permutation matrix.

Crosstalk errors given in (6) and (23) are equal while C is a generalized permutation matrix (maximum values in every row share different columns). However, the crosstalk error increases when the global separation matrix is in the situations described in (Zhang et al., 2011). The more maximum values in different rows share the same column, the worse the effect of blind source separation is and the greater the crosstalk error calculated by equation (23) is. Thus, the proposed adaptive calculation method of the crosstalk error can evaluate the effectiveness of blind source separation algorithm. Plugging the solutions of the location optimizations into the algorithm, we get the complete adaptive algorithm for calculating the crosstalk error. The pseudo code of the proposed algorithm is given in Algorithm 1.

Algorithm 1

Pseudo code of the proposed algorithm

4 Experiments and Analysis

In this section, the proposed algorithm is evaluated by two methods, which are simulation experiment and real world experiment.

4.1 Simulation Experiments

The new proposed crosstalk error definition defined by formula (23) is used to recalculate the crosstalk error of the matrixes described by the formulas (7). The new crosstalk error values are

PI (C_{1}) = 1.3000

$\mathit{PI}({\mathbf{C}_{1}})=1.3000$

and

PI (C_{2}) = 1.3021

$\mathit{PI}({\mathbf{C}_{2}})=1.3021$

, respectively. It can be seen that the new crosstalk error that can reflect the effectiveness of BSS is poor, and improves the instability of the crosstalk error that due to the small value changes effectively.

Fig. 1

Comparison of the estimation of

PI

$\mathit{PI}$

and the theoretical value of

PI

$\mathit{PI}$

In order to verify the robustness of the new proposed calculation method for evaluating BSS, two simulation examples for comparison are presented. First, we use the natural gradient RLS algorithm mentioned in literature (Zhu and Zhang, 2002) as the BSS algorithm. The source signal vector is

\begin{matrix} S_{t} = {[sign (cos (310 π t)), sin (1600 π t), sin (600 π t + 5 cos (60 π t)), sin (180 π t)]}^{T}, \end{matrix}

\[ {\mathbf{S}_{t}}={\big[\operatorname{sign}(\cos (310\pi t)),\hspace{2.5pt}\sin (1600\pi t),\hspace{2.5pt}\sin (600\pi t+5\cos (60\pi t)),\hspace{2.5pt}\sin (180\pi t)\big]^{T}},\]

where the superscript T is the meaning of matrix transposition. In this experiment, the sampling rate is taken as 10 KHz. The number of sampling points is 5000. A is randomly generated by the function

A = randn(4,4)

$\mathbf{A}=\texttt{randn(4,4)}$

of MATLAB library. Let

W = I

$\mathbf{W}=\mathbf{I}$

, where I is the unit matrix. Zero-mean independent Gaussian white noise is added to the mixtures, and the signal-to-noise ratio is 20 dB. The forgetting factor of the algorithm is 0.993. The proximity of the estimation of

PI

$\mathit{PI}$

calculated by our proposed algorithm and the theoretical value of

PI

$\mathit{PI}$

calculated by the conventional RLS algorithm is shown in Fig. 1, the value of the crosstalk error is the average of 500 Monte Carlo trials.

According to the initial conditions, the initial values of the estimated and theoretical values are 0.8669 and 0.8412, respectively. From Fig. 1, it is found that the difference between the estimated value and the theoretical value in the initial 500 iterations is obvious, but then the difference between them decreases with the increase of the number of iterations. Obviously, the trend of them is consistent and the estimated value of our proposed algorithm is more accurate than the theoretical value.

Fig. 2

The source signals.

Fig. 3

The estimation of source signals.

Figure 2 shows the source signal waveforms. Figure 3 shows the estimated signal waveforms of the source signals during the process. It can be clearly observed that the estimation of the source signals are quite similar to the source signals, which indicates the effectiveness of our proposed algorithm.

The former experiment was established under positive definite conditions. In this case, we evaluate the algorithm in over-determined conditions, and the source signal vector becomes:

\begin{matrix} S_{t} = {[sin (1600 π t), sin (180 π t)]}^{T} . \end{matrix}

\[ {\mathbf{S}_{t}}={\big[\sin (1600\pi t),\hspace{2.5pt}\sin (180\pi t)\big]^{T}}.\]

The mixed matrix becomes

A = randn(4,2)

$\mathbf{A}=\texttt{randn(4,2)}$

, and the number of samples is 16,000. The remaining experimental conditions remain unchanged. The 500 Monte Carlo trials are run, and convergence curves of the estimated value and the theoretical value of are shown in Fig. 4. As can be seen from Fig. 4, compared with the conventional algorithm, our proposed algorithm has a faster convergence rate and a smaller convergence value. Obviously, our proposed algorithm can still work well in the case of an over-determined model.

Fig. 4

Average performance index in an over-determined model.

Fig. 5

Average performance index of the proposed algorithm and the conventional algorithm.

In order to indicate the advantage of our proposed algorithm in terms of convergence rate, we apply to the conventional algorithm proposed by Zhang et al. (2011) to the same data set. The 500 Monte Carlo trials are run, and the comparison curves of crosstalk error in iterative process is shown in Fig. 5. From Fig. 5, we can see that the conventional algorithm proposed by Zhang requires about 15,000 iterations to achieve convergence, but our proposed algorithm only requires about 1800 iterations. Obviously, our proposed algorithm increases the convergence rate more than eight times relative to the conventional algorithm proposed by Zhang.

4.2 Real World Signals Experiments

In the real word signals experiments, a set of actual signals are used to test the proposed algorithm. The experimental environment is shown in Fig. 6. The signal transmission system consists of a signal generator and four transmitting antennas. The signal receiving system is made up of a signal receiver and four-element antenna array.

Fig. 6

The experiment setup.

In the test, the signal generator generates single-frequency signals that are transmitted through the transmitting antennas, and then the signal receiver receives the mixed signals using the four-element antenna array. The sampling frequency is set to 62 MHz. We select 16,000 sampling points to compute the crosstalk error. The experimental result is shown in Fig. 7. It can be clearly observed that the convergence value of the proposed algorithm is smaller than that of the conventional algorithm and the convergence rate of the proposed algorithm is far better than the conventional algorithm. It demonstrates that the proposed algorithm is effective in practical experiment.

Fig. 7

Average performance index of the proposed algorithm and the conventional algorithm in a real world experiment.

From the above experiments, we can know that the proposed algorithm for calculating the crosstalk error not only improves the robustness and the convergence rate of the crosstalk error greatly, but also deduces the crosstalk error to practical applications.

5 Conclusion

At present, the crosstalk error is widely used to verify the validity and stability of the BSS algorithms in the simulation conditions. However, in the real application environment, it is impossible to know the global separation matrix of the signal in advance. Therefore, this criterion cannot be used in practice. In this paper, the crosstalk error calculation model is established based on the RLS algorithm, and the calculation method of the crosstalk error is deduced as an adaptive algorithm without prior knowledge by using real-time estimating of the global separation matrix. It can greatly extend the use scope of the crosstalk error. At the same time, a new crosstalk error definition is proposed to improve the robustness and convergence rate of the original crosstalk error. Finally, the experimental results show that the method proposed in this paper can predict the convergence trend of the crosstalk error well, which indicates the validity and robustness of the method proposed in this paper.

Footnotes

Acknowledgements

The authors would like to thank the editor and the referees for their helpful comments and suggestions that have improved the presentation.

References

Amari

S.I.

(1997). Neural learning in structured parameter spaces-natural Riemannian gradient. In: Advances in Neural Information Processing Systems, pp. 127–133.

Aroudi

Mirkovic

De Vos

Doclo

(2016). Auditory attention decoding with EEG recordings using noisy acoustic reference signals. In: 2016 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, pp. 694–698.

Basak

Amari

S.I.

(1999). Blind separation of uniformly distributed signals: a general approach. IEEE Transactions on Neural Networks, 10(5), 1173–1185.

Bridwell

D.A.

Rachakonda

Silva

R.F.

Pearlson

G.D.

Calhoun

V.D.

(2018). Spatiospectral decomposition of multi-subject EEG: evaluating blind source separation algorithms on real and realistic simulated data. Brain Topography, 31(1), 47–61.

Choi

Cichocki

(1997). August. Adaptive blind separation of speech signals: Cocktail party problem. In: International Conference on Speech Processing, pp. 617–622.

Cichocki

Amari

S.I.

(2002). Adaptive Blind Signal and Image Processing: Learning Algorithms and Applications, Vol. 1. John Wiley and Sons.

El-Sankary

Asai

Motomura

Kuroda

(2015). Crosstalk Rejection in 3-D-Stacked Interchip Communication With Blind Source Separation. IEEE Transactions on Circuits and Systems II: Express Briefs, 62(8), 726–730.

Grellier

Comon

(1998). Blind separation of discrete sources. IEEE Signal Processing Letters, 5(8), 212–214.

Jiang

(2014). Blind signal separation for fault analysis of 400 Hz solid-state power supply. In: 2014 33rd Chinese Control Conference (CCC). IEEE, pp. 3213–3218.

10.

Laheld

Cardoso

J.F.

(1993). Separation adaptative de sources en aveugle. Implantation complexe sans contraintes. In: Actes du XIVeme Colloque GRETSI, pp. 329–332.

11.

Sejnowski

T.J.

(1995). Adaptive separation of mixed broadband sound sources with delays by a beamforming Herault–Jutten network. IEEE Journal of Oceanic Engineering, 20(1), 73–79.

12.

Feng

Z.P.

Chu

F.L.

(2012). A performance index of blind source separation based on improved crosstalk index. Zhendong yu Chongji (Journal of Vibration and Shock), 31(18), 83–88.

13.

Macchi

(1993). Self-adaptive source separation by direct and recursive networks. In: Proceedings International Conference on Digital Signal Processing (DSP’93) Limasol, Cyprus, pp. 1154–1159.

14.

Macchi

Moreau

(1999). Adaptive unsupervised separation of discrete sources. Signal Processing, 73(1–2), 49–66.

15.

Mansour

Kawamoto

Ohnishi

(2002). A survey of the performance indexes of ICA algorithms. In: Proc. IASTED Int. Conf. on Modelling, Identification and Control (MIC), pp. 660–666.

16.

Mirzal

(2017). NMF versus ICA for blind source separation. Advances in Data Analysis and Classification, 11(1), 25–48.

17.

Moreau

Macchi

(1996). High-order contrasts for self adaptive source separation. International Journal of Adaptive Control and Signal Processing, 10(1), 19–46.

18.

Moreau

Macchi

(1998). Self-adaptive source separation. II. Comparison of the direct, feedback, and mixed linear network. IEEE Transactions on Signal Processing, 46(1), 39–50.

19.

Nandakumar

M.M.

Bijoy

K.E.

(2014). October. Performance evaluation of single channel speech separation using non-negative matrix factorization. In: 2014 National Conference on Communication, Signal Processing and Networking (NCCSN). IEEE, pp. 1–4.

20.

Parra

Sajda

(2003). Blind source separation via generalized eigenvalue decomposition. Journal of Machine Learning Research, 4, 1261–1269.

21.

Qian

Chunsheng

Jinfang

(2006). Adaptive algorithm for blind source separation from instantaneous noisy mixtures. Journal of Data Acquisition and Processing, 1–8.

22.

Yang

H.H.

(1999). Serial updating rule for blind separation derived from the method of scoring. IEEE Transactions on Signal Processing, 47(8), 2279–2285.

23.

Yang

H.H.

Amari

S.I.

(1997). Adaptive online learning algorithms for blind separation: maximum entropy and minimum mutual information. Neural Computation, 9(7), 1457–1482.

24.

Yatabe

Kitamura

(2018). April. Determined blind source separation via proximal splitting algorithm. In: 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, pp. 776–780.

25.

Zhang

Gao

Hou

(2011). Use estimation of performance index to realize adaptive blind source separation. In: 2011 4th International Congress on Image and Signal Processing (CISP), Vol. 5. IEEE, pp. 2322–2326.

26.

Zhu

X.L.

Zhang

X.D.

(2002). Adaptive RLS algorithm for blind source separation using a natural gradient. IEEE Signal Processing Letters, 9(12), 432–435.