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

A kind of complex field blind source separation method

Technical field

The invention belongs to blind signal processing technology field, relate to a kind of complex field blind source separation method.

Background technology

Blind source separating (BlindSourceSeparation, BSS) also known as Blind Signal Separation (BlindSignalSeparation, BSS), be widely used in radio communication, radar, image, voice, biomedicine, the fields such as seismic wave detection are signal transacting field hot research problems.

In order to realize blind source separating, when the prioris such as source signal and aliasing parameter thereof are all unknown, fuzzy with under the fuzzy acceptable prerequisite of yardstick in arrangement, a lot of method often utilizes the aliasing signal that receives and based on the statistical property of source signal, construct the objective matrix that a group has diagonalizable structure, by carrying out Joint diagonalization operation to objective matrix group, be called as the estimation of the aliasing matrix (or its inverse matrix) of " Joint diagonalization device ", and and then Restorer varieties signal accordingly, realize blind source separating.

But existing Joint diagonalization algorithm, there is its limitation mostly, be embodied in the many restrictions to constructed objective matrix.Such as, document ^[1]require that objective matrix is positive definite matrix.Document ^[2]-[3]requirement has at least an objective matrix to be positive definite matrix, to construct whitening matrix accordingly, solves after " Joint diagonalization device " prewhitening asked is unitary matrice again.But, because objective matrix group is all obtained by statistical method usually, be limited to the impact of sample number and noise, itself there is certain error in the diagonalizable structure of each objective matrix, and the Joint diagonalization precision that pre-whitening operation is equivalent to sacrifice all the other objective matrixs has exchanged the strict diagonalization of positive definite objective matrix for, introduces error.And this error introduced in the prewhitening stage, cannot be corrected in orthogonal Joint diagonalization algorithm subsequently, thus be have impact on the overall performance of algorithm.

In view of the error that the above-mentioned prewhitening stage causes, many scholars propose some non-orthogonal joint diagonalization algorithms not needing whitening pretreatment successively, as J-Di algorithm ^[4], QDIAG algorithm ^[5], WEDGE algorithm ^[6], ACDC algorithm ^[7], SVDJD algorithm ^[8], SeDJoCo algorithm ^[9], FAJD algorithm ^[10], CVFFDIAG algorithm ^[11]deng.In the algorithm of these excellent performances, a lot of algorithm all hypothetical target matrix group is necessary for real-valued, this means that these real number field Joint diagonalization algorithms can only solve the blind source separating problem that aliasing matrix and source signal are real number value, when being applied to some key areas such as communication when blind source separating, many signals of communication are complex valued signals; In addition, when solving convolution mixed blind source separation problem, usually use easier frequency domain method, and the object to be separated on each Frequency point, often complex value, real number field Joint diagonalization algorithm cannot solve this type of complex field problem.For the algorithm that can realize complex field Joint diagonalization, ACDC algorithm and SeDJoCo algorithm require that objective matrix is hermitian symmetrical matrix or real symmetric matrix, and SVDJD algorithm requires that objective matrix is positive definite Hermitian matrix.Obviously, to the restriction of objective matrix, greatly constrain the range of application of above-mentioned algorithm.

List of references:

[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.

Summary of the invention

For above-mentioned problems of the prior art or defect, the object of the invention is to, a kind of applicability improving Joint diagonalization algorithm is provided, effectively solve the complex field blind source separation method of complex field blind source separating problem.

To achieve these goals, the present invention adopts following technical scheme:

A kind of complex field blind source separation method, comprises the following steps:

Step one: according to T observation signal sample of observation signal x (t) structure complex field objective matrix group { C ^k, k=1 ... K}, and

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

Wherein, D ^k, k=1 ... K is diagonal matrix, A ∈ F ^{m × N}(M>=N) is aliasing matrix, A ^hfor the associate matrix of aliasing matrix A;

Step 2: the complex field objective matrix group { C that step one is constructed ^k, k=1 ... K} carries out real symmetrization, the objective matrix group that the real-valued objective matrix obtaining re-constructing is formed

Step 3: utilize objective matrix group build Joint diagonalization least square cost function;

Step 4: the Joint diagonalization least square cost function utilizing alternately least-squares iteration Algorithm for Solving step 3 to obtain, obtains the estimation of aliasing matrix A according to obtain the estimation of source signal realize blind source separating, wherein represent inverse matrix, x (t) represents observation signal.

Particularly, the process of described step 2 comprises:

According to the objective matrix group C that step one constructs ^k2K matrix group is re-constructed according to formula (3) and (4) with

${\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 ^kHrepresent C ^kassociate matrix, Re (D ^k) represent diagonal matrix D ^kreal part, Im (D ^k) represent diagonal matrix D ^kimaginary part;

Note matrixing function is 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)$

By formula (3) and formula (5),

${\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 ^t(A) be the transposition of f (A), $\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];$

By formula (4) and (5),

${\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];$

By formula (6) and formula (7), the objective matrix group that the real-valued objective matrix that must re-construct is formed

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

Wherein, represent transposition; having can Joint diagonalization structure, for 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}.$

Particularly, the process of described step 3 comprises:

Build Joint diagonalization least square 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, the estimated matrix of H representing matrix A, Λ ^k, k=1 ..., 2K represents diagonal matrix group estimated matrix; represent Frobenius norm; F ^tfor estimated matrix, represent transposed matrix, and H=EF, E are Generalized Permutation Matrix.

Particularly, the process of described step 4 comprises:

Step 4.1: note iterative 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 m time iteration, Joint diagonalization least square cost function is: $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 Λ ^km () asks smallest point, obtain diagonal matrix Λ ^k(m), k=1 ..., the estimated value of 2K, wherein, H (m-1) represents the estimated value of the H that the m-1 time iteration obtains, and F (m-1) represents the estimated value of the F that the m-1 time iteration obtains;

Step 4.3: the diagonal matrix Λ that the estimated value F (m-1) of the F utilizing the m-1 time iteration to obtain and step 4.2 obtain ^k(m), k=1 ..., 2K, asks and makes Joint diagonalization least square cost function $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$ Reach the H (m) during minimum value;

Step 4.4: the Λ that the H (m) utilizing step 4.3 to obtain and step 4.2 obtain ^k(m), k=1 ..., 2K, asks and makes $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$ Reach the F (m) during minimum value;

Step 4.5: whether evaluation algorithm restrains, namely judges whether to meet and wherein δ is the threshold value of setting; If do not met, return step 4.2; If met, then obtain the estimation of final H and F, the two all can be used as the estimation of matrix A, according to obtain the estimation of aliasing matrix A according to obtain the estimation of source signal realize blind source separating, wherein represent inverse matrix, x (t) represents observation signal.

Compared with prior art, the present invention has following technique effect:

1, complex field objective matrix group is constructed, and carry out real symmetrization, the objective matrix group that the real-valued objective matrix obtaining re-constructing is formed, is converted into real number field Joint diagonalization problem by complex field Joint diagonalization problem, may be used for solving complex field blind source separating problem; And compared with other are equally applicable to complex field algorithm, do not retrain and treat that Joint diagonalization objective matrix is hermitian symmetrical matrix or positive definite Hermitian matrix, there is extremely wide applicability.

2, adopt the alternately least-squares iteration algorithm based on Joint diagonalization least square cost function, make full use of the design feature of the objective matrix group that real-valued objective matrix is formed, to realize the Joint diagonalization of fresh target matrix group.

3, utilize alternately least-squares iteration Algorithm for Solving cost function, obtain the estimated value of aliasing matrix, realize blind source separating, simulation results show the inventive method has higher convergence precision.

4, present invention eliminates pre-whitening operation, avoid and introduce albefaction error.

Accompanying drawing explanation

Fig. 1 is the flow chart of the inventive method;

Fig. 2 is that in experiment one, the overall situation is refused to make an uproar the performance curve that horizontal GRL changes with iterations;

Fig. 3 is the performance curve that in experiment one, parameter JH changes with iterations

Fig. 4 is under signal to noise ratio NER different in experiment two, and GRL is with the change curve of iterations;

Fig. 5 is under signal to noise ratio NER different in experiment two, and average GRL is with the change curve of iterations;

Fig. 6 adopts distinct methods, and when signal to noise ratio NER is 15dB, average GRL is with the change curve of iterations;

Fig. 7 adopts distinct methods, and when signal to noise ratio NER is 10dB, average GRL is with the change curve of iterations;

Fig. 8 is the change curve of average GRL with iterations of experiment 3 100 independent experiments;

Fig. 9 is the planisphere of 4 source signals;

Figure 10 is the planisphere of 4 Received signal strength;

Figure 11 is the planisphere of 4 restoring signals;

Below in conjunction with drawings and Examples, explanation detailed further and explanation are done to the solution of the present invention.

Embodiment

Defer to technique scheme, complex field blind source separation method of the present invention, see Fig. 1, specifically comprises the following steps:

Step one: structure objective matrix group.

The linear array be made up of M array element, N number of narrow-band source incides in linear array from different directions, and the M that linear array receives ties up observation signal x (t) and is:

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

Wherein, A ∈ F ^{m × N}(M>=N) is aliasing matrix, s (t)=[s ₁(t) ..., s _n(t)] ^tfor source signal vector, x (t)=[x ₁(t) ..., x _m(t)] ^tfor observation signal, n (t)=[n ₁(t) ..., n _m(t)] ^tfor noise vector.

According to T observation signal sample of observation signal x (t) construct the objective matrix group { C that one group of total number is K ^k, k=1 ... K}, wherein each objective matrix all has following diagonalizable structure:

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

Wherein D ^k, k=1 ... K is diagonal matrix, A ^hrepresent the associate matrix of aliasing matrix A.

Consider the restriction of sample number and the existence of noise, formula (2) described diagonal structure is only approximate existence.The main method of structure objective matrix group has: the covariance matrix under time windows ^[1], the cross-correlation matrix under different time shift ^[2], the section matrix under three rank or more Higher Order Cumulants ^[12], spatial time-frequency matrix ^[13], or second feature function ^[14].Method of the present invention can be used for solving the most general complex field Joint diagonalization problem to realize complex field blind source separating, does not have particular/special requirement to objective matrix, namely allows aliasing matrix A and diagonal matrix D ^k, k=1 ... K is complex valued matrices, and therefore the building method of objective matrix group can select the combination of above-mentioned a kind of or any several method.

Step 2: real symmetrization is carried out to the objective matrix group that step one constructs, the objective matrix group that the real-valued objective matrix obtaining re-constructing is formed complex field Joint diagonalization problem is converted into real number field Joint diagonalization problem.

According to the objective matrix group C that step one constructs ^k2K matrix group is re-constructed according to formula (3) and (4) with

${\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 ^kHrepresent C ^kassociate matrix, Re (D ^k) represent diagonal matrix D ^kreal part, Im (D ^k) represent diagonal matrix D ^kimaginary part.

Note matrixing function is 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)$

By formula (3) and formula (5),

${\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 ^t(A) be the transposition of f (A), $\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].$

By formula (4) and formula (5),

${\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].$

By formula (6) and formula (7), the objective matrix group that the real-valued objective matrix that must re-construct is formed

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

K dimension is the matrix group { C that the complex value objective matrix of M × M is formed by above-mentioned conversion process ^k, k=1 ... K} is converted into the objective matrix group that 2K dimension is the real-valued objective matrix formation that 2M × 2M is new find by analyzing, new objective matrix group there is following design feature:

(T1) objective matrix group having can Joint diagonalization structure, namely for diagonal matrix, even row vector ${\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 [] represents that structure take vector as the diagonal matrix of diagonal entry.

(T2) objective matrix is real symmetric matrix, namely

(T3), in expression formula (8), the representation 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)$

From formula (9), if obtained front M column element, utilize its design feature, can be easy to obtain m column element below.

Step 3: the objective matrix group utilizing step 2 to obtain build Joint diagonalization least square cost function.

Complex field general goals matrix group is converted into real number field symmetric targets matrix group by the matrixing process described in step 2, for formula (8), can utilize the real number field Joint diagonalization algorithm of existing excellent performance, recovers matrix A, realizes blind source separating.But, because existing real number field Joint diagonalization algorithm is unspecial in problem formula (8) Suo Shi, thus cannot take into full account parametric structures feature described in (T1) ~ (T3) and be used in the optimizing process of algorithm.In order to make full use of these design features as priori in optimizing process to improve the performance of algorithm, the present invention proposes the alternately least-squares iteration algorithm based on Joint diagonalization least square cost function, to realize the Joint diagonalization of fresh target matrix group.

The cost function form of conventional sign Joint diagonalization degree mainly contains: information theory criterion cost function, F norm minimum criterion cost function and least square cost function.According to formula (8), the Joint diagonalization least square cost function that the present invention adopts, is shown below:

$\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 representing matrix estimated matrix, H ^trepresenting matrix transposed matrix estimated matrix, Λ ^k, k=1 ..., 2K represents diagonal matrix group estimated matrix, if diagonal matrix group with represent, then ${\tilde{d}}^k=\[{\widetilde{\mathrm{\λ}}}_1^k,\mathrm{...},{\widetilde{\mathrm{\λ}}}_M^k\],{\widetilde{\mathrm{\λ}}}_i^k,i=1,\mathrm{...},M$ Represent described in step 2 (T1) estimated value; represent Frobenius norm;

In formula (10), objective matrix group be known, ask for H and Λ ^k, k=1 ..., the estimated value of 2K, makes H Λ ^kh ^t, k=1 ..., 2K approaches as much as possible that is, cost function is made little as much as possible, be the object of following steps of the process.By observing cost function known, it is biquadratic function for unknown matrix H, and thus, computation complexity is higher, not easily asks for function smallest point.

For the ease of calculating, cost function is handled 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 above formula, introduce matrix F ^tfor estimated matrix, cost function is deteriorated to respectively about H and F by the biquadratic function about unknown matrix H ^tquadratic function, reduce and solve difficulty.And, utilize known H and F of the character of matrix decomposition be essence equal, i.e. H=EF, wherein E is Generalized Permutation Matrix, the every a line namely in E each show and only have a nonzero element.And, H and F ^tobviously (T3) described matrix in rapid two should all be had design feature, this design feature is used in follow-up optimizing process.

Step 4: utilize alternately least-squares iteration Algorithm for Solving cost function.

Introduce a kind of alternately least-squares iteration algorithm based on gradient descent method, this alternative and iterative algorithm each take turns iteration and can be divided into three steps, in each step, alternately estimate H and F, and diagonal matrix Λ ^k, k=1 ..., 2K, the smallest point of search cost function J, and the design feature making full use of involved each parameter in the process.

Step 4.1: according to the design feature of (T3) described matrix A in step 2, note iterative 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 m time iteration, Joint diagonalization least square cost function is: $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 Λ ^km () asks smallest point, obtain diagonal matrix Λ ^k(m), k=1 ..., the estimated value of 2K, wherein, H (m-1) represents the estimated value of the H that the m-1 time iteration obtains, and F (m-1) represents the estimated value of the F that the m-1 time iteration obtains.Concrete methods of realizing is as follows:

To 2K Joint diagonalization least square cost function, respectively about Λ ^kthe individual different diagonal entry of M of (m) differentiate, and make derivative be zero, obtain row vector described in step 2 (T1) estimation

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

Wherein, (i, the j) of matrix P individual element p _ijfor:

$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 ^ki-th element be:

$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 formula (13) and (14), h _ii-th column vector of representing matrix H (m-1), parameter h _i+M, h _j, h _j+Mcan reasoning accordingly; f _ii-th column vector of representing matrix F (m-1), f _i, f _i+M, f _j, f _j+Mimplication can reasoning accordingly.So far,

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

Step 4.3: the diagonal matrix Λ that the estimated value F (m-1) of the F utilizing the m-1 time iteration to obtain and step 4.2 obtain ^k(m), k=1 ..., 2K, asks and makes Joint diagonalization least square cost function $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$ Reach the H (m) during minimum value, concrete methods of realizing is as follows:

Note H (m)=[H ₁(m), H ₂(m)], wherein H ₁(m) and H ₂m the 1st row that () is respectively H (m) arrange 2M column vector to M column vector and M+1, 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, with be respectively G ^km 1 to the M row vector of () and M+1 are to 2M row vector;

To Joint diagonalization least square cost function about H ₁m () differentiate also makes derivative be zero, due to the m time iteration H ₂m () not yet upgrades, still use the iteration result H that m-1 takes turns ₂(m-1), 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 by H ₁m (), according to (T3) described matrix in step 2 $\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}$ Design feature, can H be obtained ₂(m), and then obtain H (m).

Step 4.4: the Λ that the H (m) utilizing step 4.3 to obtain and step 4.2 obtain ^k(m), k=1 ..., 2K, asks and makes $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$ Reach the F (m) during minimum value, concrete methods of realizing is as follows:

Note F (m)=[F ₁(m), F ₂(m)], wherein F ₁(m) and F ₂m the 1st row that () is respectively F (m) arrange 2M column vector to M column vector and M+1, and wherein with be respectively L ^km 1 to the M column vector of () and M+1 are to 2M column vector;

To Joint diagonalization least square cost function about F ₁m () differentiate also makes derivative be zero, in computational process, due to F ₂m () not yet upgrades, still use the iteration result F of the m-1 time ₂(m-1), 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 (), according to (T3) described matrix in step 2 $\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}$ Design feature, can F be obtained ₂(m), and then obtain F (m).

Step 4.5: whether evaluation algorithm restrains, namely judges whether to meet and wherein δ is the threshold value of setting, is less positive number, is set to 0.01; If not, step 4.2 is returned; If so, then obtain the estimation of final H and F, the two difference is very little, all can be used as matrix estimation, according to obtain the estimation of aliasing matrix A according to obtain the estimation of source signal realize blind source separating, wherein represent inverse matrix, x (t) represents observation signal.

Emulation experiment:

Define two performance index, first performance digit synbol is JH, is 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)$

Above-mentioned parameter JH is obviously consistent with the cost function value constructed by the inventive method, the variation tendency of the inventive method cost function value in an iterative process of reflection, because parameter JH fails the estimated performance that directly embodies aliasing matrix, therefore be only applied in experiment one, so that the validity of the inventive method for the cost function shown in formula (11) to be described.

Second individual character energy index is called that the overall situation refuses the level of making an uproar (GlobalRejectionLevel, GRL), in order to weigh aliasing Matrix Estimation value and the difference between actual value A, be 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 _ijrepresent matrix the i-th row jth column element.

Experiment one: the validity and the convergence that utilize not Noise objective matrix group verification algorithm.

Run 20 independent experiments, each experiment all constructs the not Noise objective matrix group { C of N × N ^k=A Λ ^ka ^h, k=1 ... K}, wherein K=15, N=10, aliasing matrix A and diagonal matrix Λ ^kbe the random complex valued matrices produced.The overall situation that what Fig. 2 provided is is refused to make an uproar the performance curve that horizontal GRL changes with iterations, what abscissa represented is iterations, what ordinate represented is that the overall situation refuses the horizontal GRL that makes an uproar, in order to obtain this curve, need in each iteration (taking turns iteration for m), the relational expression according to H (m) convolution (9) of trying to achieve $\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},$ Obtain the aliasing Matrix Estimation value of current iteration and then the GRL of current iteration is tried to achieve according to formula (18).What Fig. 3 provided is the performance curve that parameter JH changes with iterations.Fig. 2 and Fig. 3 all illustrates, when noiseless and objective matrix have accurately can Joint diagonalization structure, algorithm of the present invention reaches convergence with high precision.As shown in Figure 1, after algorithmic statement, aliasing Matrix Estimation value with difference between actual value is minimum, algorithm of the present invention effectively can realize the estimation of aliasing matrix.

Experiment two: utilize the band of structure to make an uproar the rapidity of objective matrix group verification algorithm and validity.

Provide the target square formation group of N × N:

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

Wherein K=15, N=10.Aliasing matrix A and diagonal matrix Λ ^kbe the random complex valued matrices produced.In order to characterize the intensity of noise matrix, incite somebody to action not Noise part A Λ ^ka ^hwith noise section Δ C ^kratio be 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 simulations, random generation complex values noise matrix Δ C ^k(k=1 ..., K) and to meet NER=10dB respectively, 15dB, 20dB, 25dB.Under different signal to noise ratio NER values, run 20 independent experiments respectively, each independently realize GRL with iterations change curve as shown in Figure 4, under different N ER value, average GRL is shown in Fig. 5 with the change curve of iterations, two figure all show under situation made an uproar by band, and namely objective matrix is that under non-critical diagonalizable structure situation, algorithm of carrying still can be restrained with degree of precision.

Adopt CVFFDIAG method respectively ^[11], FAJD method ^[10]15dB is respectively at signal to noise ratio NER with the inventive method (STBJD), when 20dB, respectively run 20 independent experiments, average GRL with iterations change curve as shown in Figure 6 and Figure 7, compared with the inventive method and all the other two kinds of methods, can restrain with degree of precision.

Experiment three: the validity of checking the inventive method when solving blind source separating problem.

Provide the complex value source signal of 4 zero-mean statistical iteration: 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), herein suppose that 4 information sources are received by 4 array elements, aliasing channel parameter is the random complex valued matrices A produced, and sample of signal number T=400, the noise introduced in receiving course is N=randn (4, T)+jrandn (4, T).Signal-to-noise ratio settings is 10dB.Objective matrix is obtained and selected number K=8 by the cross-correlation matrix asked under the different time shift of Received signal strength.Run the estimation of algorithm realization aliasing matrix A of the present invention, and Restorer varieties signal.

The average GRL of 100 independent experiments with iterations change curve as shown in Figure 8.In 100 independent experiments, randomly draw the 33rd experiment, check the recovery effects of source signal: Fig. 9-Figure 11 represents the planisphere of 4 source signals respectively, the planisphere of 4 Received signal strength, and the planisphere of 4 restoring signals.As seen from the figure, source signal is fully mixed by aliasing matrix and under there is larger noise situations, method of the present invention still can realize convergence and convergence error is less, achieves the comparatively accurate blind separation of source signal.Meanwhile, restoring signal shown in source signal and Figure 11 shown in comparison diagram 9, clearly can find the arrangement ambiguity that blind source separating field is intrinsic and yardstick ambiguity phenomenon.

Claims (4)

1. a complex field blind source separation method, is characterized in that, comprises the following steps:

Step one: according to T observation signal sample of observation signal x (t) structure complex field objective matrix group { C ^k, k=1 ... K}, and

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

Wherein, D ^k, k=1 ... K is diagonal matrix, A ∈ F ^{m × N}(M>=N) is aliasing matrix, A ^hfor the associate matrix of aliasing matrix A;

Step 2: the complex field objective matrix group { C that step one is constructed ^k, k=1 ... K} carries out real symmetrization, the objective matrix group that the real-valued objective matrix obtaining re-constructing is formed

Step 3: utilize objective matrix group build Joint diagonalization least square cost function;

Step 4: the Joint diagonalization least square cost function utilizing alternately least-squares iteration Algorithm for Solving step 3 to obtain, obtains the estimation of aliasing matrix A according to obtain the estimation of source signal realize blind source separating, wherein represent inverse matrix, x (t) represents observation signal.

2. complex field blind source separation method as claimed in claim 1, it is characterized in that, the process of described step 2 comprises:

According to the objective matrix group C that step one constructs ^k2K matrix group is re-constructed according to formula (3) and (4) with

${\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 ^kHrepresent C ^kassociate matrix, Re (D ^k) represent diagonal matrix D ^kreal part, Im (D ^k) represent diagonal matrix D ^kimaginary part;

Note matrixing function is 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)$

By formula (3) and formula (5),

${\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 ^t(A) be the transposition of f (A), $\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];$

By formula (4) and (5),

${\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];$

By formula (6) and formula (7), the objective matrix group that the real-valued objective matrix that must re-construct is formed

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

Wherein, represent transposition; having can Joint diagonalization structure, k=1 ..., 2K is 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. complex field blind source separation method as claimed in claim 2, it is characterized in that, the process of described step 3 comprises:

Build Joint diagonalization least square 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 representing matrix estimated matrix, Λ ^k, k=1 ..., 2K represents diagonal matrix group k=1 ..., the estimated matrix of 2K; represent Frobenius norm; F ^tfor estimated matrix, represent transposed matrix, and H=EF, E are Generalized Permutation Matrix.

4. complex field blind source separation method as claimed in claim 3, it is characterized in that, the process of described step 4 comprises:

Step 4.1: note iterative 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 m time iteration, Joint diagonalization least square cost function is: $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,$ K=1 ..., 2K, about Λ ^km () asks smallest point, obtain diagonal matrix Λ ^k(m), k=1 ..., the estimated value of 2K, wherein, H (m-1) represents the estimated value of the H that the m-1 time iteration obtains, and F (m-1) represents the estimated value of the F that the m-1 time iteration obtains;

Step 4.3: the diagonal matrix Λ that the estimated value F (m-1) of the F utilizing the m-1 time iteration to obtain and step 4.2 obtain ^k(m), k=1 ..., 2K, asks and makes Joint diagonalization least square cost function $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$ Reach the H (m) during minimum value;

Step 4.4: the Λ that the H (m) utilizing step 4.3 to obtain and step 4.2 obtain ^k(m), k=1 ..., 2K, asks and makes $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$ Reach the F (m) during minimum value;

Step 4.5: whether evaluation algorithm restrains, namely judges whether to meet and wherein δ is the threshold value of setting; If do not met, return step 4.2; If met, then obtain the estimation of final H and F, the two all can be used as matrix estimation, according to obtain the estimation of aliasing matrix A according to obtain the estimation of source signal realize blind source separating, wherein represent inverse matrix, x (t) represents observation signal.