CN105282067A - Complex field blind source separation method - Google Patents

Abstract

The invention discloses a complex field blind source separation method. A complex filed target matrix system is built, and real symmetrization is carried out to obtain a reconstructed target matrix system formed by a real-value target matrix, the complex field combined diagonalization problem is converted into the real field combined diagonalization problem to solve the complex field blind source separation problem; compared with other algorithms that are also suitable for the complex filed, the method doesn't restrain a diagonalization target matrix to be combined into a hermitian symmetric matrix or a positive definite hermitian matrix and is wide in application; an alternative least square iterative algorithm based on combined diagonalization least square cost functions is adopted, and the structural characteristics of the target matrix system formed by the real-value target matrix are fully used to realize the combined diagonalization of a new target matrix system; The cost functions are solved by the alternative least square iterative algorithm, the estimate values of a mixed matrix are obtained, the blind source separation is realized, and the simulation results verify that the method provided is high in convergence precision.

Description

Complex field blind source separation method

Technical Field

The invention belongs to the technical field of blind signal processing, and relates to a complex field blind source separation method.

Background

The Blind Source Separation (BSS), also known as Blind Signal Separation (BSS), is widely applied to the fields of wireless communication, radar, image, voice, biomedicine, seismic wave detection and the like, and is a hot research topic in the field of signal processing.

In order to realize blind source separation, under the condition that prior knowledge of a source signal, aliasing parameters thereof and the like is unknown, and under the premise that permutation ambiguity and scale ambiguity are acceptable, a plurality of methods often utilize the received aliasing signal and based on the statistical characteristics of the source signal to construct a group of target matrixes with diagonalizable structures, and perform joint diagonalization operation on the group of target matrixes to obtain the estimation of the aliasing matrix (or the inverse matrix thereof) called as a joint diagonalizer, and further recover the source signal according to the estimation, so as to realize blind source separation.

However, most of the existing joint diagonalization algorithms have limitations, particularly in terms of limitations on the constructed target matrix. For example, the literature^[1]The target matrices are required to be positive definite matrices. Literature reference^[2]-[3]At least one target matrix is required to be a positive definite matrix so as to construct a whitening matrix, and a joint diagonalizer to be solved is pre-whitened into a unitary matrix and then solved. However, since the target matrix set is usually obtained by a statistical method, and is limited by the number of samples and the influence of noise, the diagonalizable structure of each target matrix has a certain error, and the pre-whitening operation is equivalent to strict diagonalization of a positive target matrix at the expense of the joint diagonalization accuracy of the rest of target matrices, and an error is introduced. And the error introduced in the pre-whitening stage cannot be corrected in the subsequent orthogonal joint diagonalization algorithm, thereby affecting the overall performance of the algorithm.

In view of the error caused by the pre-whitening stage, many scholars have proposed several non-orthogonal joint diagonalization algorithms, such as the J-Di algorithm, without whitening pre-processing^[4]QDIAG algorithm^[5]Wedge algorithm^[6]ACDC algorithm^[7]SVDJD algorithm^[8]SeDJoco algorithm^[9]The FAJD algorithm^[10]CVFFDIAG algorithm^[11]And the like. In the algorithms with excellent performance, a target matrix group is assumed to be real-valued by a plurality of algorithms, which means that the real-domain joint diagonalization algorithms can only solve the blind source separation problem that the aliasing matrix and the source signals are both real-valued, and when the blind source separation is applied to some important fields such as communication, a plurality of communication signals are complex-valued signals; in addition, when solving the problem of the separation of the convolution mixed blind source, a simpler frequency domain method is often used, and the object to be separated on each frequency point is often complex-valued, and the complex-valued domain joint diagonalization algorithm cannot solve the problem of the complex number domain. For the algorithm capable of realizing complex domain joint diagonalization, the ACDC algorithm and the SeDJoco algorithm require that the target matrix is a hermitian symmetric matrix or a real symmetric matrix, and the SVDJD algorithm requires that the target matrix is a positive definite hermitian matrix. Obviously, the limitation on the target matrix greatly restricts the application range of the algorithm.

Reference documents:

[1]D.-T.Pham.Jointapproximatediagonalizationofpositivedefinitematrices[J].SIAMJ.MatrixAnal.Appl.,2001,55(4):1136–1152.

[2]A.Belouchrani，K.A.Meraim，J.-F.Cardoso，etal.Ablindsourceseparationtechniqueusingsecond-orderstatistics[J].IEEETrans.SignalProcess.,1997,45(2):434-444.

[3]M.Wax,J.Sheinvald.Aleast-squaresapproachtojointdiagonalization[J].IEEESignalProcess.Lett.,1997,4(2):52–53.

[4]A.Souloumiac.Nonorthogonaljointdiagonalizationbycombininggivensandhyperbolicrotations[J].IEEETrans.SignalProcess.,2009,57(6):2222-2231.

[5]R.Vollgraf,K.Obermayer.Quadraticoptimizationforsimultaneousmatrixdiagonalization[J].IEEETrans.SignalProcess.,2006,54(9):3270-3278.

[6]P.A.Yeredor.Fastapproximatejointdiagonalizationincorporatingweightmatrices[J].IEEETrans.SignalProcess.,2009,57(3):878-891.

[7]A.Yeredor.Non-orthogonaljointdiagonalizationintheleast-squaressensewithapplicationinblindsourceseparation[J].IEEETrans.SignalProcess.,2002,50(7):1545-1553.

[8]K.Todros,J.Tabrikian.QML-BasedJointDiagonalizationofPositive-DefiniteHermitianMatrices*J+.IEEETrans.SignalProcess.,2010,58(9):4656-4673.

[9]A.Yeredor,B.Song,F.Roemer,etal.Asequentiallydrilledjointcongruence(SeDJoCo)transformationwithapplicationsinblindsourceseparationandmultiuserMIMOsystems[J].IEEETrans.SignalProcess.,2012,60(6):2744-2757.

[10]X.-L.Li,X.-D.Zhang.Nonorthogonaljointdiagonalizationfreeofdegeneratesolution[J].IEEETrans.SignalProcess.,2007,55(5):1803-1814.

[11]X.-F.Xu,D.-Z.Feng,W.X.Zheng.Afastalgorithmfornonunitaryjointdiagonalizationanditsapplicationtoblindsourceseparation[J].IEEETrans.SignalProcess.,2011,59(7):3457-3463.

[12]X.-L.Li,X.-D.Zhang.Nonorthogonaljointdiagonalizationfreeofdegeneratesolution[J].IEEETrans.SignalProcess.,2007,55(5):1803-1814.

[13]X.-F.Xu,D.-Z.Feng,W.X.Zheng.Afastalgorithmfornonunitaryjointdiagonalizationanditsapplicationtoblindsourceseparation[J].IEEETrans.SignalProcess.,2011,59(7):3457-3463.

[14]P.Zhang,S.Chen,L.Hanzo.Embeddediterativesemi-blindchannelestimationforthree-stage-concatenatedMIMO-aidedQAMturbotransceivers[J].IEEETrans.Veh.Technol.,2014,61(1):439–446.

disclosure of Invention

In view of the above problems or drawbacks of the prior art, an object of the present invention is to provide a complex field blind source separation method that can improve the applicability of the joint diagonalization algorithm and effectively solve the problem of complex field blind source separation.

In order to achieve the purpose, the invention adopts the following technical scheme:

a complex field blind source separation method comprises the following steps:

the method comprises the following steps: t observation signal samples from the observation signal x (T)Constructing a complex field object matrix set { C^kK is 1, … K }, and

C^k＝AD^kA^H,k＝1,…K(2)

wherein D is^kK1, … K is a diagonal matrix, a ∈ F^M×N(M.gtoreq.N) is a aliasing matrix, A^HA conjugate transpose matrix of the aliasing matrix A;

step two: for the complex field object matrix set { C constructed in step one^kK is 1, … K, and the real symmetry is carried out to obtain an object matrix group composed of reconstructed real value object matrix

Step three: using sets of objective matricesConstructing a joint diagonalization least square cost function;

step four: solving the joint diagonalization least square cost function obtained in the step three by using an alternative least square iterative algorithm to obtain the estimation of the aliasing matrix AAccording toObtaining an estimate of a source signalRealize blind source separation, whereinTo representX (t) represents the observed signal.

Specifically, the process of the second step includes:

constructing a target matrix set C according to the step one^kReconstructing 2K matrix groups according to equations (3) and (4)And

${\tilde{C}}^{2k}=\frac{1}{2}(C^k+C^{kH})=A\mathrm{Re}\left(D^k\right)A^H,k=1,...,K---\left(3\right)$

${\tilde{C}}^{2k-1}=\frac{1}{2j}(C^k-C^{kH})=A\mathrm{Im}\left(D^k\right)A^H,k=1,...,K---\left(4\right)$

wherein, C^kHIs represented by C^kConjugate transpose matrix of Re (D)^k) Representing a diagonal matrix D^kReal part of, Im (D)^k) Representing a diagonal matrix D^kAn imaginary part of (d);

let the matrix transformation function be f (X):

$f\left(X\right)=\left[\begin{array}{cc}\mathrm{Re}\left(X\right)& \mathrm{Im}\left(X\right)\\ -\mathrm{Im}\left(X\right)& \mathrm{Re}\left(X\right)\end{array}\right]---\left(5\right)$

from formula (3) and formula (5) to obtain

${\stackrel{\‾}{C}}^{2k}\stackrel{\mathrm{\Δ}}{=}f\left({\tilde{C}}^{2k}\right)=f\left(A\right)\left[\begin{array}{cc}\mathrm{Re}\left(D^k\right)& 0\\ 0& \mathrm{Re}\left(D^k\right)\end{array}\right]f^T\left(A\right)=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^{2k}{\stackrel{\‾}{A}}^T,k=1,...,K---\left(6\right)$

Wherein f is^T(A) As a result of the transpose of (f), $\stackrel{\‾}{A}=f\left(A\right),{\stackrel{\‾}{A}}^T=f\left(A\right),{\stackrel{\‾}{D}}^{2k}=\left[\begin{array}{cc}\mathrm{Re}\left(D^k\right)& 0\\ 0& \mathrm{Re}\left(D^k\right)\end{array}\right];$

from the formulae (4) and (5) to give

${\stackrel{\‾}{C}}^{2k-1}\stackrel{\mathrm{\Δ}}{=}f\left({\tilde{C}}^{2k-1}\right)=f\left(A\right)\left[\begin{array}{cc}\mathrm{Im}\left(D^k\right)& 0\\ 0& \mathrm{Im}\left(D^k\right)\end{array}\right]f^T\left(A\right)=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^{2k-1}{\stackrel{\‾}{A}}^T,k=1,\mathrm{...},K---\left(7\right)$

Wherein, ${\stackrel{\‾}{D}}^{2k-1}=\left[\begin{array}{cc}\mathrm{Im}\left(D^k\right)& 0\\ 0& \mathrm{Im}\left(D^k\right)\end{array}\right];$

target matrix group composed of reconstructed real-valued target matrix from equations (6) and (7)

${\stackrel{\‾}{C}}^k=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^k{\stackrel{\‾}{A}}^T,k=1,...,2K---\left(8\right)$

Wherein,to representTransposing;has a structure which can be combined with diagonalization,is a diagonal matrix; $\stackrel{\‾}{A}=f\left(A\right)=\left[\begin{array}{cc}\mathrm{Re}\left(A\right)& \mathrm{Im}\left(A\right)\\ -\mathrm{Im}\left(A\right)& \mathrm{Re}\left(A\right)\end{array}\right]\∈R^{2M\×2M}.$

specifically, the process of the third step includes:

construction of a Joint diagonalized least squares cost function J (H, Λ)¹,…,Λ^2K,F)：

$J(H,{\mathrm{\Λ}}^1,...,{\mathrm{\Λ}}^{2K},F)=\underset{k=1}{\overset{2K}{\Σ}}\left|\right|{\stackrel{\‾}{C}}^k-{\mathrm{H\Λ}}^k{F}^T|{|}_F^2---\left(11\right)$

Where H represents the estimated matrix of matrix A, Λ^kK1, …,2K denotes a diagonal array groupThe estimation matrix of (2);represents the Frobenius norm; f^TIs composed ofIs determined by the estimation matrix of (a),to representAnd H ═ EF, E is the generalized permutation matrix.

Specifically, the process of the fourth step includes:

step 4.1, recording iteration initial value $H\left(0\right)=F\left(0\right)=\left[\begin{array}{cc}{I}_{M\×M}& I_{M\×M}\\ -I_{M\×M}& I_{M\×M}\end{array}\right];$

Step 4.2, in the mth iteration, the joint diagonalization least square cost function is as follows: $J_k\left(H\right(m-1),{\mathrm{\Λ}}^k(m),F(m-1\left)\right)=\left|\right|{\stackrel{\‾}{C}}^k-H(m-1){\mathrm{\Λ}}^k\left(m\right)F^T(m-1)|{|}_F^2,k=1,...,2K,$ about Λ^k(m) minimum to obtain diagonal array Λ^k(m), K is an estimated value of 1, …,2K, where H (m-1) denotes the (m-1) th iterationF (m-1) represents the estimated value of F obtained in the (m-1) th iteration;

step 4.3, utilizing the estimated value F (m-1) of F obtained in the (m-1) th iteration and the diagonal array Λ obtained in the step 4.2^k(m), K1, …,2K, such that a joint diagonalized least squares cost function is solved $J\left(H\right(m),{\mathrm{\Λ}}^1(m),\mathrm{...},{\mathrm{\Λ}}^{2K}(m),F(m-1\left)\right)=\underset{k=1}{\overset{2K}{\Σ}}\left|\right|{\stackrel{\‾}{C}}^k-H\left(m\right){\mathrm{\Λ}}^k\left(m\right)F^T(m-1)|{|}_F^2$ H (m) when the minimum value is reached;

step 4.4: utilizing the steps4.3H (m) and Λ from step 4.2^k(m), K is 1, …,2K, in the formula $J\left(H\right(m),{\mathrm{\Λ}}^1(m),...,{\mathrm{\Λ}}^{2K}\left(m\right),F\left(m\right))=\underset{k=1}{\overset{2K}{\Σ}}||{\stackrel{\‾}{C}}^k-H\left(m\right){\mathrm{\Λ}}^k\left(m\right)F^T\left(m\right)|{|}_F^2$ F (m) when the minimum value is reached;

step 4.5: judging whether the algorithm converges, i.e. whether it satisfiesAnd isWherein is a set threshold value; if it is notIf not, returning to the step 4.2; if so, a final estimate of H and F is obtained, both of which can be used as estimates of the matrix A, based on Obtaining an estimate of an aliasing matrix AAccording toObtaining an estimate of a source signalRealize blind source separation, whereinTo representX (t) represents the observed signal.

Compared with the prior art, the invention has the following technical effects:

1. constructing a complex field target matrix group, carrying out real symmetry to obtain a target matrix group formed by a reconstructed real value target matrix, converting the complex field joint diagonalization problem into a real field joint diagonalization problem, and solving the complex field blind source separation problem; compared with other algorithms which are also suitable for complex field, the target matrix to be united and diagonalized is not restricted to be Hermite symmetric matrix or positive definite Hermite matrix, and the method has extremely wide applicability.

2. And an alternating least square iterative algorithm based on a joint diagonalization least square cost function is adopted, and the structural characteristics of a target matrix group formed by real-valued target matrices are fully utilized to realize the joint diagonalization of the new target matrix group.

3. The cost function is solved by using the alternating least square iterative algorithm to obtain the estimated value of the aliasing matrix, the blind source separation is realized, and the simulation result verifies that the method has higher convergence precision.

4. The invention cancels the pre-whitening operation and avoids the introduction of whitening error.

Drawings

FIG. 1 is a flow chart of the method of the present invention;

FIG. 2 is a performance curve of global noise rejection level GRL as a function of iteration number in experiment one;

FIG. 3 is a performance curve of parameter JH varying with iteration number in experiment I

FIG. 4 is a graph of GRL versus iteration number for different signal-to-noise ratios NER in experiment two;

FIG. 5 is a graph of the average GRL versus iteration number for different signal-to-noise ratios NER in experiment two;

FIG. 6 is a graph of the average GRL versus iteration number using a different approach, with a signal-to-noise ratio NER of 15 dB;

FIG. 7 is a graph of the average GRL versus iteration number using a different approach, with a signal-to-noise ratio NER of 10 dB;

FIG. 8 is a plot of average GRL versus iteration number for three 100 independent experiments;

fig. 9 is a constellation diagram of 4 source signals;

fig. 10 is a constellation of 4 received signals;

fig. 11 is a constellation diagram of 4 recovered signals;

the invention will be explained and explained in more detail below with reference to the drawings and exemplary embodiments.

Detailed Description

According to the above technical solution, the complex field blind source separation method of the present invention, referring to fig. 1, specifically includes the following steps:

the method comprises the following steps: an objective matrix set is constructed.

The array comprises M array elements, N narrow-band sources are incident on the array from different directions, and M-dimensional observation signals x (t) received by the array are as follows:

x(t)＝As(t)+n(t)(1)

wherein, A ∈ F^M×N(M.gtoreq.N) is an aliasing matrix, s (t) [ s ]₁(t),…,s_N(t)]^TIs a vector of source signals, x (t) ═ x₁(t),…,x_M(t)]^TTo observe the signal, n (t) ═ n₁(t),…,n_M(t)]^TIs a noise vector.

T observation signal samples from the observation signal x (T)Constructing a set of target matrix set { C with total number K^kK1, … K, where each object matrix has a diagonalizable structure as follows:

C^k＝AD^kA^H,k＝1,…K(2)

wherein D^kK is 1, … K is a diagonal matrix, a^HRepresenting the conjugate transpose of the aliasing matrix a.

Considering the limitation of the number of samples and the existence of noise, the diagonalization structure in equation (2) is only approximately present. The main methods for constructing the target matrix set are: covariance matrices under different time windows^[1]Cross correlation matrix at different time shifts^[2]Slice matrix at third or higher order cumulants^[12]Space time-frequency matrix^[13]Or second characteristicFunction(s)^[14]. The method of the invention can be used for solving the most common problem of complex field joint diagonalization to realize complex field blind source separation, has no special requirement on a target matrix, namely allows a aliasing matrix A and a diagonal matrix D^kAnd K is 1, … K is a complex value matrix, so that the construction method of the target matrix group can select one or a combination of any two of the methods.

Step two: the target matrix group constructed in the step one is subjected to real symmetry to obtain a target matrix group constructed by a reconstructed real-value target matrixThe complex domain joint diagonalization problem is converted into a real domain joint diagonalization problem.

Constructing a target matrix set C according to the step one^kReconstructing 2K matrix groups according to equations (3) and (4)And

${\tilde{C}}^{2k}=\frac{1}{2}(C^k+C^{kH})=A\mathrm{Re}\left(D^k\right)A^H,k=1,...,K---\left(3\right)$

${\tilde{C}}^{2k-1}=\frac{1}{2j}(C^k-C^{kH})=A\mathrm{Im}\left(D^k\right)A^H,k=1,...,K---\left(4\right)$

wherein, C^kHIs represented by C^kConjugate transpose matrix of Re (D)^k) Representing a diagonal matrix D^kReal part of, Im (D)^k) Representing a diagonal matrix D^kThe imaginary part of (c).

Let the matrix transformation function be f (X):

$f\left(X\right)=\left[\begin{array}{cc}\mathrm{Re}\left(X\right)& \mathrm{Im}\left(X\right)\\ -\mathrm{Im}\left(X\right)& \mathrm{Re}\left(X\right)\end{array}\right]---\left(5\right)$

from formula (3) and formula (5) to obtain

${\stackrel{\‾}{C}}^{2k}\stackrel{\mathrm{\Δ}}{=}f\left({\tilde{C}}^{2k}\right)=f\left(A\right)\left[\begin{array}{cc}\mathrm{Re}\left(D^k\right)& 0\\ 0& \mathrm{Re}\left(D^k\right)\end{array}\right]f^T\left(A\right)=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^{2k}{\stackrel{\‾}{A}}^T,k=1,...,K---\left(6\right)$

Wherein f is^T(A) As a result of the transpose of (f), $\stackrel{\‾}{A}=f\left(A\right),{\stackrel{\‾}{A}}^T=f\left(A\right),{\stackrel{\‾}{D}}^{2k}=\left[\begin{array}{cc}\mathrm{Re}\left(D^k\right)& 0\\ 0& \mathrm{Re}\left(D^k\right)\end{array}\right].$

from formula (4) and formula (5) to obtain

${\stackrel{\‾}{C}}^{2k-1}\stackrel{\mathrm{\Δ}}{=}f\left({\tilde{C}}^{2k-1}\right)=f\left(A\right)\left[\begin{array}{cc}\mathrm{Im}\left(D^k\right)& 0\\ 0& \mathrm{Im}\left(D^k\right)\end{array}\right]f^T\left(A\right)=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^{2k-1}{\stackrel{\‾}{A}}^T,k=1,\mathrm{...},K---\left(7\right)$

Wherein, ${\stackrel{\‾}{D}}^{2k-1}=\left[\begin{array}{cc}\mathrm{Im}\left(D^k\right)& 0\\ 0& \mathrm{Im}\left(D^k\right)\end{array}\right].$

target matrix group composed of reconstructed real-valued target matrix from equations (6) and (7)

${\stackrel{\‾}{C}}^k=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^k{\stackrel{\‾}{A}}^T,k=1,...,2K.---\left(8\right)$

The transformation process described above uses a set of K complex-valued object matrices of dimension M × M { C }^kK1, … K into a target matrix set of 2K new real-valued target matrices of dimension 2M × 2MThrough analysis finding, a new target matrix setHas the following structural characteristics:

(T1) target matrix setWith structures which can be jointly diagonalized, i.e.As diagonal matrices, i.e. if row vectors ${\stackrel{\‾}{d}}^k=\[{\mathrm{\λ}}_1^k,\mathrm{...},{\mathrm{\λ}}_M^k\]\∈R^{1\×M},$ Then ${\stackrel{\‾}{D}}^k=diag\[{\stackrel{\‾}{d}}^k,{\stackrel{\‾}{d}}^k\]\∈R^{2M\×2M},$ Wherein diag [. C]The expression constructs a diagonal matrix with vector · as diagonal elements.

(T2) the target matrices are all real symmetric matrices, i.e.

(T3) in expression (8), the expression of matrix a is as follows:

$\stackrel{\‾}{A}=f\left(A\right)=\left[\begin{array}{cc}\mathrm{Re}\left(A\right)& \mathrm{Im}\left(A\right)\\ -\mathrm{Im}\left(A\right)& \mathrm{Re}\left(A\right)\end{array}\right]\∈R^{2M\×2M}---\left(9\right)$

as is clear from formula (9), if obtainedThe first M rows of elements can be easily obtained by utilizing the structural characteristicsThe latter M columns of elements.

Step three: utilizing the target matrix group obtained in the step twoAnd constructing a joint diagonalized least square cost function.

And the matrix transformation process converts the complex number domain common target matrix group into a real number domain symmetrical target matrix group, and can recover the matrix A by using the existing real number domain joint diagonalization algorithm with excellent performance according to the formula (8) to realize blind source separation. However, since the existing real-number domain joint diagonalization algorithm does not specifically address the problem shown in equation (8), the parametric structural features described in (T1) to (T3) cannot be fully considered and utilized in the optimization process of the algorithm. In order to fully utilize the structural characteristics as prior knowledge in the optimization process to improve the performance of the algorithm, the invention provides an alternating least square iterative algorithm based on a joint diagonalization least square cost function to realize the joint diagonalization of a new target matrix group.

The commonly used form of the cost function characterizing the degree of joint diagonalization is mainly: an information theory criterion cost function, an F-norm minimization criterion cost function, and a least squares cost function. According to equation (8), the joint diagonalized least squares cost function employed by the present invention is shown as follows:

$\stackrel{\‾}{J}(H,{\mathrm{\Λ}}^1,...,{\mathrm{\Λ}}^{2K})=\underset{k=1}{\overset{2K}{\Σ}}\left|\right|{\stackrel{\‾}{C}}^k-{\mathrm{H\Λ}}^k{H}^T|{|}_F^2---\left(10\right)$

wherein H represents a matrixOf the estimation matrix, H^TRepresentation matrixIs transposed matrix ofΛ of the estimation matrix^kK1, …,2K denotes a diagonal array groupIf the diagonal array group isBy usingIs shown to be ${\tilde{d}}^k=\[{\widetilde{\mathrm{\λ}}}_1^k,\mathrm{...},{\widetilde{\mathrm{\λ}}}_M^k\],{\widetilde{\mathrm{\λ}}}_i^k,i=1,\mathrm{...},M$ Denotes that described in step two (T1)An estimated value of (d);represents the Frobenius norm;

in the formula (10), the object matrix groupIt is known to find H and Λ^kK is an estimate of 1, …,2K, such that H Λ^kH^TK is as close as possible to 1, …,2KI.e. make the cost functionAs small as possible, is the object of the following steps of the invention. By observing a cost functionIt can be known that it is a quartic function for the unknown matrix H, and therefore, the computational complexity is high, and it is not easy to find the minimum point of the function.

For the convenience of calculation, the cost function is processed as follows:

$J(H,{\mathrm{\Λ}}^1,\mathrm{...},{\mathrm{\Λ}}^{2K},F)=\underset{k=1}{\overset{2K}{\Σ}}\left|\right|{\stackrel{\‾}{C}}^k-{\mathrm{H\Λ}}^k{F}^T|{|}_F^2---\left(11\right)$

in the above formula, a matrix F is introduced^TIs composed ofSuch that the cost function is degenerated from a quartic function with respect to the unknown matrix H to a quartic function with respect to H and F, respectively^TThe quadratic function of (2) reduces the solving difficulty. Furthermore, by the nature of matrix decomposition, H and F are essentially equal, i.e., H ═ EF, where E is a generalized permutation matrix, i.e., each row in E has one and only one non-zero element per column. And, H and F^TIt is obvious that each should have the matrix described in step two (T3)The structural characteristics of (a) will be utilized in the subsequent optimization process.

Step four: and solving the cost function by using an alternating least square iterative algorithm.

An alternative least square iterative algorithm based on gradient descent method is introduced, which can be divided into three steps in each iteration, H and F are estimated alternatively in each step, and a diagonal matrix Λ is formed^kAnd K is 1, … and 2K, searching the minimum point of the cost function J, and fully utilizing the structural characteristics of the parameters involved in the process.

Step 4.1, recording an iteration initial value according to the structural characteristics of the matrix A in the step two (T3) $H\left(0\right)=F\left(0\right)=\left[\begin{array}{cc}{I}_{M\×M}& I_{M\×M}\\ -I_{M\×M}& I_{M\×M}\end{array}\right];$

Step 4.2, in the mth iteration, the joint diagonalization least square cost function is as follows: $J_k\left(H\right(m-1),{\mathrm{\Λ}}^k(m),F(m-1\left)\right)=\left|\right|{\stackrel{\‾}{C}}^k-H(m-1){\mathrm{\Λ}}^k\left(m\right)F^T(m-1)|{|}_F^2,k=1,...,2K,$ about Λ^k(m) minimum to obtain diagonal array Λ^kAnd (m), K is an estimated value of 1, … and 2K, wherein H (m-1) represents an estimated value of H obtained in the (m-1) th iteration, and F (m-1) represents an estimated value of F obtained in the (m-1) th iteration. The specific implementation method comprises the following steps:

for 2K joint diagonalized least squares cost functions, Λ is applied^k(M) M distinct diagonal elementsTaking the derivative and making the derivative zero to obtain the line vector in step two (T1)Is estimated by

${\tilde{d}}^k={\left(P^{-1}{q}^k\right)}^T---\left(12\right)$

Wherein the (i, j) th element P of the matrix P_ijComprises the following steps:

$p_{ij}=h_i^T{h}_j{f}_j^T{f}_i+h_i^T{h}_{j+M}{f}_{j+M}^T{f}_i+h_{i+M}^T{h}_j{f}_j^T{f}_{i+M}+h_{i+M}^T{h}_{j+M}{f}_{j+M}^T{f}_{i+M};---\left(13\right)$

column vector q^kThe ith element of (a) is:

$q_i^k=h_i^T{\stackrel{\‾}{C}}^k{f}_i+h_{i+M}^T{\stackrel{\‾}{C}}^k{f}_{i+M}.---\left(14\right)$

in formulae (13) and (14), h_iThe ith column vector, parameter H, representing the matrix H (m-1)_i+M,h_j,h_j+MCan reason accordingly; f. of_iThe ith column vector, F, representing the matrix F (m-1)_i,f_i+M,f_j,f_j+MThe meaning of (c) can be inferred accordingly. At this point in time,

then diagonal matrix ${\mathrm{\Λ}}^k\left(m\right)=diag\[{\tilde{d}}^k,{\tilde{d}}^k\].$

Step 4.3, utilizing the estimated value F (m-1) of F obtained in the (m-1) th iteration and the diagonal array Λ obtained in the step 4.2^k(m), K1, …,2K, such that a joint diagonalized least squares cost function is solved $J\left(H\right(m),{\mathrm{\Λ}}^1(m),\mathrm{...},{\mathrm{\Λ}}^{2K}(m),F(m-1\left)\right)=\underset{k=1}{\overset{2K}{\Σ}}\left|\right|{\stackrel{\‾}{C}}^k-H\left(m\right){\mathrm{\Λ}}^k\left(m\right)F^T(m-1)|{|}_F^2$ H (m) when the minimum value is reached, the specific implementation method is as follows:

h (m) ═ H₁(m),H₂(m)]In which H is₁(m) and H₂(M) column 1 to M column vectors and M +1 to 2M column vectors of H (M), respectively, and ${\mathrm{\Λ}}^k\left(m\right)F^T(m-1)=G^k\left(m\right)=\left[\begin{array}{c}{G}_1^k\left(m\right)\\ G_2^k\left(m\right)\end{array}\right],$ wherein,andare each G^k(M) the 1 st to mth row vectors and the M +1 st to 2 mth row vectors;

for joint diagonalization least square cost function with respect to H₁(m) taking the derivative and making the derivative zero, since the mth iteration H₂(m) not yet updatedStill using the iteration result H of round m-1₂(m-1) to obtain:

$H_1\left(m\right)=\[\left(\underset{k=1}{\overset{2K}{\Σ}}{\stackrel{\‾}{C}}^k{G}_1^{kT}\left(m\right)\right)-H_2\left(m-1\right)\left(\underset{k=1}{\overset{2K}{\Σ}}{G}_2^k\right(m\left)G_1^{kT}\left(m\right)\right)\]{\left(\underset{k=1}{\overset{2K}{\Σ}}{G}_1^k\right(m\left)G_1^{kT}\left(m\right)\right)}^{-1}---\left(15\right)$

then is formed by H₁(m) said matrix according to step two (T3) $\stackrel{\‾}{A}=f\left(A\right)=\left[\begin{array}{cc}\mathrm{Re}\left(A\right)& \mathrm{Im}\left(A\right)\\ -\mathrm{Im}\left(A\right)& \mathrm{Re}\left(A\right)\end{array}\right]\∈R^{2M\×2M}$ Structural feature of (1), H can be obtained₂(m) to obtain H (m).

Step 4.4 utilization of H (m) from step 4.3 and Λ from step 4.2^k(m), K is 1, …,2K, in the formula $J\left(H\right(m),{\mathrm{\Λ}}^1(m),\mathrm{...},{\mathrm{\Λ}}^{2K}(m),F(m\left)\right)=\underset{k=1}{\overset{2K}{\Σ}}\left|\right|{\stackrel{\‾}{C}}^k-H\left(m\right){\mathrm{\Λ}}^k\left(m\right)F^T(m-1)|{|}_F^2$ F (m) when the minimum value is reached, the specific implementation method is as follows:

notation F (m) ═ F₁(m),F₂(m)]In which F is₁(m) and F₂(M) column 1 to M column vectors and M +1 to 2M column vectors, respectively, of F (M), andwhereinAndand are each L^k(M) 1 st to M th column vectors and M +1 st to 2M th column vectors;

for joint diagonalization least squares cost function with respect to F₁(m) taking the derivative and making the derivative zero, during the calculation, due to F₂(m) not yet updated, still using the m-1 st iteration result F₂(m-1) to obtain:

$F_1\left(m\right)=\[\left(\underset{k=1}{\overset{2K}{\Σ}}{\stackrel{\‾}{C}}^{kT}{L}_1^k\left(m\right)\right)-F_2\left(m-1\right)\left(\underset{k=1}{\overset{2K}{\Σ}}{L}_2^{kT}\right(m\left)L_1^k\left(m\right)\right)\]{\left(\underset{k=1}{\overset{2K}{\Σ}}{L}_1^{kT}\right(m\left)L_1^k\left(m\right)\right)}^{-1}---\left(16\right)$

then by F₁(m) said matrix according to step two (T3) $\stackrel{\‾}{A}=f\left(A\right)=\left[\begin{array}{cc}\mathrm{Re}\left(A\right)& \mathrm{Im}\left(A\right)\\ -\mathrm{Im}\left(A\right)& \mathrm{Re}\left(A\right)\end{array}\right]\∈R^{2M\×2M}$ Structural feature of (1), F can be obtained₂(m) to obtain F (m).

Step 4.5: judging whether the algorithm converges, i.e. whether it satisfiesAnd isWherein the threshold value is set to be a smaller positive number and is set to be 0.01; if not, returning to the step 4.2; if so, then the final H andthe estimate of F, with little difference, can be used as a matrixIs estimated according toObtaining an estimate of an aliasing matrix AAccording toObtaining an estimate of a source signalRealize blind source separation, whereinTo representX (t) represents the observed signal.

Simulation experiment:

two performance indicators are defined, the first one being JH, defined as:

$JH=10\lg \[\underset{k=1}{\overset{2K}{\Σ}}\left|\right|{\stackrel{\‾}{C}}^k-{\mathrm{H\Λ}}^k{F}^T|{|}_F^2\]---\left(17\right)$

the parameter JH is obviously consistent with the cost function value constructed by the method of the present invention, and reflects the variation trend of the cost function value in the iterative process of the method of the present invention, and since the parameter JH cannot directly reflect the estimation performance of the aliasing matrix, the parameter JH is only applied to the first experiment to illustrate the effectiveness of the method of the present invention for the cost function shown in formula (11).

The second performance measure is called global noise rejection level (GRL), which is used to measure the aliasing matrix estimation valueAnd the true value a, defined as:

$GRL=\underset{i=1}{\overset{N}{\Σ}}\left(\underset{j=1}{\overset{N}{\Σ}}\frac{\left|g_{ij}\right|}{{\max }_k\left|g_{ik}\right|}-1\right)+\underset{j=1}{\overset{N}{\Σ}}\left(\underset{i=1}{\overset{N}{\Σ}}\frac{\left|g_{ij}\right|}{{\max }_k\left|g_{kj}\right|}-1\right)---\left(18\right)$

wherein, g_ijRepresentative matrixRow i and column j.

Experiment one: and verifying the effectiveness and convergence of the algorithm by using the noise-free target matrix set.

20 independent experiments were run, each to construct a set of N × N noise-free target matrices { C }^k＝AΛ^kA^HK1, … K, where K15 and N10, aliasing matrix a and diagonal matrix Λ^kAre all randomly generated complex value matrices. FIG. 2 shows a global noise rejection level GRL as a function ofIn order to obtain a performance curve with a variable iteration number, the abscissa represents the iteration number, and the ordinate represents the global noise rejection level GRL, it is necessary to obtain a performance curve with a variable iteration number according to the relation shown in the obtained h (m) and the combination formula (9) in each iteration (in the case of the mth iteration) $\stackrel{\‾}{A}=f\left(A\right)=\left[\begin{array}{cc}\mathrm{Re}\left(A\right)& \mathrm{Im}\left(A\right)\\ -\mathrm{Im}\left(A\right)& \mathrm{Re}\left(A\right)\end{array}\right]\∈R^{2M\×2M},$ Obtaining the aliasing matrix estimation value of the iterationFurther, the GRL of the current iteration is obtained from equation (18). FIG. 3 shows the parameter JH as a function of the number of iterationsPerformance curves. Both fig. 2 and fig. 3 illustrate that the algorithm of the present invention converges with a very high accuracy in the case of no noise, i.e. the target matrix has an accurately joint diagonalizing structure. As can be seen from FIG. 1, after the algorithm converges, the aliasing matrix estimatesThe difference between the real value and the real value is extremely small, and the algorithm of the invention can effectively realize the estimation of the aliasing matrix.

Experiment two: and verifying the rapidity and the effectiveness of the algorithm by using the constructed noisy target matrix group.

Given the target square array set of nxn:

C^k＝AΛ^kA^H+ΔC^k(k＝1,…,K)(19)

where K15, N10, aliasing matrix a and diagonal matrix Λ^kAll are randomly generated complex value matrices, to characterize the strength of the noise matrix, the noise free part a Λ^kA^HAnd noise part Δ C^kThe ratio of (d) is expressed as NER:

$NER=10{\log }_{10}\frac{\left|\right|{\mathrm{A\Λ}}^k{A}^H|{|}_F^2}{\left|\right|{\mathrm{\ΔC}}^k|{|}_F^2},k=1,...,K---\left(20\right)$

in the simulation, a complex-valued noise matrix Δ C is randomly generated^k(K1, …, K) to satisfy NER 10dB, 15dB, 20dB, 25dB, respectively. Under different signal-to-noise ratio NER values, 20 independent experiments are respectively run, the change curve of GRL along with the iteration times is independently realized each time as shown in fig. 4, under different NER values, the change curve of average GRL along with the iteration times is shown in fig. 5, and both the graphs show that the algorithm can still be converged at higher precision under the condition of noise, namely under the condition that a target matrix is a non-strict diagonalizable structure.

Respectively using CVFFDIAG method^[11]FAJD method^[10]And the method (STBJD) of the invention runs 20 independent experiments respectively under the condition that the signal-to-noise ratio NER is 15dB and 20dB respectively, and the change curve of the average GRL along with the iteration times is shown in figure 6 and figure 7.

Experiment three: the effectiveness of the method in solving the blind source separation problem is verified.

4 zero-mean statistically independent complex-valued source signals are given: s₁(t)＝sin(310πt)+jcos(100πt),s₂(t)＝sin(180πt)+jsin(400πt),s₃(t)＝sin(20πt)sin(600πt)+jcos(20πt)cos(600πt),s₄(t)＝sin[600πt+6cos(120πt)]+ jcos (900 π t), hereAssuming that 4 sources are received by 4 array elements, the aliasing channel parameter is a complex value matrix a generated randomly, the number of signal samples T is 400, and the noise introduced during the receiving process is N ═ randn (4, T) + jnrandn (4, T). The signal-to-noise ratio is set to 10 dB. The target matrix is obtained by obtaining cross-correlation matrixes of the received signals under different time shifts, and the selected number K is 8. The algorithm of the invention is operated to realize the estimation of the aliasing matrix A, andthe source signal is recovered.

The average GRL over the number of iterations for 100 independent experiments is shown in fig. 8. Randomly draw 33 th experiment in 100 independent experiments, and check the recovery effect of the source signal: fig. 9-11 show the constellation of 4 source signals, the constellation of 4 received signals, and the constellation of 4 recovered signals, respectively. It can be seen from the figure that, under the condition that the source signals are fully mixed by the mixing and overlapping matrix and relatively large noise exists, the method of the invention can still realize convergence and has relatively small convergence error, and relatively accurate blind separation of the source signals is realized. Meanwhile, comparing the source signal shown in fig. 9 with the recovered signal shown in fig. 11, it is possible to clearly find the permutation uncertainty and scale uncertainty phenomenon inherent in the blind source separation field.

Claims (4)

1. A complex field blind source separation method is characterized by comprising the following steps:

the method comprises the following steps: t observation signal samples from the observation signal x (T)Constructing a complex field object matrix set { C^kK is 1.. K }, and

C^k＝AD^kA^H,k＝1,...K(2)

wherein D is^kK is the diagonal momentArray, A ∈ F^M×N(M.gtoreq.N) is a aliasing matrix, A^HA conjugate transpose matrix of the aliasing matrix A;

step two: for the complex field object matrix set { C constructed in step one^kK is 1, K is real symmetrical to obtain target matrix group formed by reconstructed real value target matrix

Step three: using sets of objective matricesConstructing a joint diagonalization least square cost function;

step four: solving the joint diagonalization least square cost function obtained in the step three by using an alternative least square iterative algorithm to obtain the estimation of the aliasing matrix AAccording toObtaining an estimate of a source signalRealize blind source separation, whereinTo representX (t) represents the observed signal.

2. The complex-domain blind source separation method as claimed in claim 1, wherein the process of step two comprises:

constructing a target matrix set C according to the step one^kReconstructing 2K matrix groups according to equations (3) and (4)And

${\tilde{C}}^{2k}=\frac{1}{2}(C^k+C^{kH})=A\mathrm{Re}\left(D^k\right)A^H,k=1,...,K---\left(3\right)$

${\tilde{C}}^{2k-1}=\frac{1}{2j}(C^k-C^{kH})=A\mathrm{Im}\left(D^k\right)A^H,k=1,...,K---\left(4\right)$

wherein, C^kHIs represented by C^kConjugate transpose matrix of Re (D)^k) Representing a diagonal matrix D^kReal part of, Im (D)^k) Representing a diagonal matrix D^kAn imaginary part of (d);

let the matrix transformation function be f (X):

$f\left(X\right)=\left[\begin{array}{cc}\mathrm{Re}\left(X\right)& \mathrm{Im}\left(X\right)\\ -\mathrm{Im}\left(X\right)& \mathrm{Re}\left(X\right)\end{array}\right]---\left(5\right)$

from formula (3) and formula (5) to obtain

${\stackrel{\‾}{C}}^{2k}\stackrel{\mathrm{\Δ}}{=}f\left({\tilde{C}}^{2k}\right)=f\left(A\right)\left[\begin{array}{cc}\mathrm{Re}\left(D^k\right)& \mathrm{\θ}\\ \mathrm{\θ}& \mathrm{Re}\left(D^k\right)\end{array}\right]f^T\left(A\right)=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^{2k}{\stackrel{\‾}{A}}^T,k=1,...,K---\left(6\right)$

Wherein f is^T(A) As a result of the transpose of (f), $\stackrel{\‾}{A}=f\left(A\right),$ ${\stackrel{\‾}{A}}^T=f\left(A\right),$ ${\stackrel{\‾}{D}}^{2k}=\left[\begin{array}{cc}\mathrm{Re}\left(D^k\right)& \mathrm{\θ}\\ \mathrm{\θ}& \mathrm{Re}\left(D^k\right)\end{array}\right];$

from the formulae (4) and (5) to give

${\stackrel{\‾}{C}}^{2k-1}\stackrel{\mathrm{\Δ}}{=}f\left({\tilde{C}}^{2k-1}\right)=f\left(A\right)\left[\begin{array}{cc}\mathrm{Im}\left(D^k\right)& 0\\ 0& \mathrm{Im}\left(D^k\right)\end{array}\right]f^T\left(A\right)=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^{2k-1}{\stackrel{\‾}{A}}^T,k=1,...,K---\left(7\right)$

Wherein, ${\stackrel{\‾}{D}}^{2k-1}=\left[\begin{array}{cc}\mathrm{Im}\left(D^k\right)& 0\\ 0& \mathrm{Im}\left(D^k\right)\end{array}\right];$

target matrix group composed of reconstructed real-valued target matrix from equations (6) and (7)

${\stackrel{\‾}{C}}^k=\stackrel{\‾}{A}{\stackrel{\‾}{D}}^k{\stackrel{\‾}{A}}^T,k=1,...,2K---\left(8\right)$

Wherein, to representTransposing;has a structure which can be combined with diagonalization,k is 1, 2K is a diagonal matrix; $\stackrel{\‾}{A}=f\left(A\right)=\left[\begin{array}{cc}\mathrm{Re}\left(A\right)& \mathrm{Im}\left(A\right)\\ -\mathrm{Im}\left(A\right)& \mathrm{Re}\left(A\right)\end{array}\right]\∈R^{2M\×2M}.$

3. the complex-domain blind source separation method as claimed in claim 2, wherein the process of step three comprises:

construction of a Joint diagonalized least squares cost function J (H, Λ)¹,...,Λ^2K,F)：

$J(H,{\mathrm{\Λ}}^1,...,{\mathrm{\Λ}}^{2K},F)=\underset{k=1}{\overset{2K}{\Σ}}\left|\right|{\stackrel{\‾}{C}}^k-{\mathrm{H\Λ}}^k{F}^T|{|}_F^2---\left(11\right)$

Wherein H represents a matrixΛ of the estimation matrix^kK1, 2K denotes a diagonal array groupAn estimation matrix of K1.., 2K;represents the Frobenius norm; f^TIs composed ofIs determined by the estimation matrix of (a),to representAnd H ═ EF, E is the generalized permutation matrix.

4. The complex-domain blind source separation method of claim 3, wherein the process of step four comprises:

step 4.1, recording iteration initial value $H\left(0\right)=F\left(0\right)=\left[\begin{array}{cc}{I}_{M\×M}& I_{M\×M}\\ -I_{M\×M}& I_{M\×M}\end{array}\right];$

Step 4.2, in the mth iteration, the joint diagonalization least square cost function is as follows: $J_k(H(m-1),{\mathrm{\Λ}}^k\left(m\right),F(m-1))={\left|\right|{\stackrel{\‾}{C}}^k-H(m-1){\mathrm{\Λ}}^k\left(m\right)F^T(m-1)\left|\right|}_F^2,$ 1, 2K for Λ^k(m) minimum to obtain diagonal array Λ^k(m), K being 1.., 2K, where H (m-1) denotes an estimate of H obtained at the m-1 st iteration and F (m-1) denotes an estimate of F obtained at the m-1 st iteration;

step 4.3, utilizing the estimated value F (m-1) of F obtained in the (m-1) th iteration and the diagonal array Λ obtained in the step 4.2^k(m), K1, 2K, solving a cost function that results in joint diagonalization of the least squares $J(H\left(m\right),{\mathrm{\Λ}}^1\left(m\right),...,{\mathrm{\Λ}}^{2K}\left(m\right),F(m-1))=\underset{k=1}{\overset{2K}{\mathrm{\Σ}}}{\left|\right|{\stackrel{\‾}{C}}^k-H\left(m\right){\mathrm{\Λ}}^k\left(m\right)F^T(m-1)\left|\right|}_F^2$ H (m) when the minimum value is reached;

step 4.4 utilization of H (m) from step 4.3 and Λ from step 4.2^k(m), K1, 2K, in a manner such that $J(H\left(m\right),{\mathrm{\Λ}}^1\left(m\right),...,{\mathrm{\Λ}}^{2K}\left(m\right),F\left(m\right))=\underset{k=1}{\overset{2K}{\mathrm{\Σ}}}{\left|\right|{\stackrel{\‾}{C}}^k-H\left(m\right){\mathrm{\Λ}}^k\left(m\right)F^T\left(m\right)\left|\right|}_F^2$ F (m) when the minimum value is reached;

step 4.5: judging whether the algorithm converges, i.e. whether it satisfiesAnd isWherein is a set threshold value; if not, returning to the step 4.2; if so, final estimates of H and F are obtained, both of which can be used as matricesIs estimated according to Obtaining an estimate of an aliasing matrix AAccording toObtaining an estimate of a source signalRealize blind source separation, whereinTo representX (t) represents the observed signal.