CN102568493A - Underdetermined blind source separation (UBSS) method based on maximum matrix diagonal rate - Google Patents

Abstract

The invention discloses an underdetermined blind source separation (UBSS) method based on a maximum matrix diagonal rate. The method comprises the following steps of: constructing inverse matrixes of C*M/N M*M-dimensional sub matrixes of a mixed matrix (wherein M and N are respectively the number of sensors and the number of source signals); multiplying the inverse matrixes by observation signal vectors to acquire initial estimation signal vectors; and sequentially calculating the covariance matrix, the solid matrix, the absolute value matrix and the diagonal rate of each initial estimation signal vector, selecting the initial estimation signal vector corresponding to the maximum diagonal rate as estimation of a source signal vector, and thus realizing underdetermined separation of source signals. By the method, the requirement for source signal sparseness is reduced, aliasing of road source signals is allowed at each time frequency point at most, and the underdetermined separation problem of music signals and noise signals is solved. The requirement for the statistical property of the source signals is low, and the underdetermined separation problem of Gaussian signals and related signals is solved. In addition, by the method, processing of each time frequency point and each sub matrix can be executed in parallel, and hardware implementation is facilitated.

Description

A kind ofly owe to decide blind separating method based on maximum matrix diagonal angle rate

Technical field

The present invention relates to a kind of method of owing under the stable condition instantaneous mixed signal is carried out blind separation.This method separable sparse, a little less than sparse or coherent signal, field such as can be applied to signal Processing, biomedicine and communicate by letter.

Background technology

(Blind Source Separation is exactly to confirm a certain conversion according to the blended data vector that observes BSS), with the technology of recovery original signal or information source to the separation of blind source.In typical case, the observation data vector is the output of one group of sensor, and what wherein each sensor received is the various combination of source signal.Term " blind " has double implication: a. source signal can not be observed; B. how to mix be unknown to source signal.When the number of source signal during more than the number of observation signal, (this kind situation more is close to practical application for underdetermined BSS, UBSS) problem, also is the technological difficulties that linear instantaneous mixes blind separation simultaneously, has challenge in order to owe blind surely separation.

Up to the present, emerged the multiple blind source separation method of owing to decide.I.F.Gorodnitsky in 1992 etc. combined electroencephalography (eeg) inverse problem proposed the focal solution of under determined system (focal underdetermined system solver, FOCUSS), this is the method for a kind of posteriority iteration weighting; 1998, S.S.Chen proposed based on the theoretical basic back tracking method of linear programming (Basis Pursuit, BP).Owing to decide blind source separation method roughly can be divided into and owe to decide change overdetermination method, time-frequency masking method and two-step approach three major types.Owe to decide the observation signal that change overdetermination method makes new advances through structure, will owe to decide separation and become the overdetermination separation, thereby reach the purpose that source signal separates.The advantage of this method is to utilize traditional overdetermination isolation technics, and shortcoming is the quality that separating effect depends on the new observation signal that constructs, in case this signal configuration is bad, disintegrate-quality will variation.Two-step approach is actually the popularization of basic back tracking method, is through separating linear equation through obtaining optimum solution to separating to use restraint, utilize the sparse property of signal, minimizes 0 norm and then retrains separating, and 0 norm deals with very inconvenient and responsive especially to noise.1999; D.L.Donoho has proved and has adopted the minimum and 0 norm minimum equivalence under certain condition of 1 norm, and 1 norm is handled well than 0 norm, utilizes linear programming can obtain optimum solution easily; Though noise robustness is better than 0 norm criterion; But effect still can not be satisfactory, and this algorithm is to be prerequisite with signal in the sparse property of time domain in addition, therefore general poor of separating.The time-frequency masking method is in 2000 the earliest; Proposed by Sam T.Roweis, 2004, Yilmaz combined the DUET algorithm to further develop the time-frequency masking algorithm with Rickard; Yet this class methods strict demand source signal is similar in whole time-frequency domain and satisfies W-dislocation orthogonality, and condition is harsh.2005, Abrard proposed the TIFROM algorithm, only had a source signal to exist in the observation signal of several time frequency windows that this method requirement is contiguous, and the length of time frequency window is difficult to confirm, and can not searches out single source time-frequency domain of all any sizes.2008, people such as M.S.Pedersen combined ICA to realize separating of voice with the time-frequency macking technique.These class methods are based on the sparse property of the time-frequency domain condition of source signal; Design a shelter template that is equivalent to Time frequency Filter; Utilize template to extract the time frequency of each source signal one by one; Thereby realized separation, the main shortcoming of this type algorithm is very strict to the requirement of the sparse property of signal, and antimierophonic poor performance.These methods all have a common defective just to be based on the sparse property of signal, and this has just limited their application greatly.In recent years; Caused scholars' concern gradually based on the weak sparse blind surely source separation problem of owing; 2007, separated in the blind surely source of owing a little less than people such as A.Aissa-El-Bey utilize the method for subspace to realize under the sparse condition, but this method operand haves much room for improvement very greatly.2009; People such as doctor Peng Dezhong utilize subspace method separation source signal under the situation that reduces the requirement of sparse property, and in 2010, have studied and do not considered the sparse property of signal; Utilize spatial time-frequency to distribute and owe the method that separate in blind surely source; The method can be separated 2*m-1 road source signal, and m is the number of instantaneous mixed signal, wherein m >=3.2011, Zhou Guo permitted to wait the people under the not strict sparse situation of signal, to propose a kind of new method-non-linear projection and row shelter that (nonlinear projection and column masking NPCM) comes estimated mixing matrix; Separate in the blind surely source of owing a little less than people such as land wind ripple have realized based on diagonalization of matrix under the sparse condition, but this method can only be handled unrelated signal.

Summary of the invention

In owing blind surely source separation problem, utilize the sparse property of source signal usually, adopt statistical method to come the separation source signal.If signal does not satisfy sparse property, or is correlated with between the signal, then separating effect is not good.To this problem, the present invention propose a kind of based on maximum matrix diagonal angle rate owe to decide blind separating method, sparse and coherent signal separated a little less than this method was applicable to.This method is intended to reduce the degree of dependence of separation algorithm to the sparse property of signal, allows signal to have aliasing to a certain degree.On the basis that the hypothesis hybrid matrix has been estimated to obtain; The present invention is at first through

individual M * M dimension submatrix of constructing hybrid matrix and the inverse matrix (wherein M, N are respectively sensor number and source signal number) of asking submatrix; Estimate signal vector at the beginning of then inverse matrix and observation signal multiplication of vectors to be obtained

individual; Last each first covariance matrix, real part matrix, absolute value matrix of estimating signal vector that calculate successively; And calculating diagonal angle rate; Get the corresponding first estimation of estimating signal vector as source signal of maximum diagonal angle rate, according to this source signal is separated.Theoretical analysis and experimental result show; Absolute value at the determinant of all submatrixs differs under the prerequisite that is not very big; Method proposed by the invention has excellent performance with respect to additive method: a. has reduced the requirement to the sparse property of source signal; Can solve the aliasing of maximum M road source signal, what solved music signal and noise signal owes to decide separation problem.B. less demanding to the statistical property of source signal, what solved gaussian signal and coherent signal owes to decide separation problem.C. the present invention is equally applicable to the very little situation of angle between the hybrid matrix column vector.Processing of frequency and each sub-matrices can executed in parallel when d. the present invention was to each, was beneficial to hardware and realized.

In order to achieve the above object, the invention provides a kind of based on maximum matrix diagonal angle rate owe to decide blind separating method, may further comprise the steps:

Step 100: utilize the observation signal of the instantaneous mixing of sensor, be expressed as:

x(t)＝As(t) (1)

X (t) in the formula (1)=[x ₁(t), x ₂(t) ..., x _M(t)] ^TThe M dimension observation signal vector that expression is received by sensor, A are represented M * N (the normalized capable full rank hybrid matrix of dimension of M＜N), s (t)=[s ₁(t), s ₂(t) ..., s _N(t)] ^TThe source signal vector that expression N dimension is unknown.

Step 200; The observation signal of instantaneous mixing is sent into the orthogonal transformation module carries out Short Time Fourier Transform, obtain time-frequency domain observation signal X (τ, w); And with this time-frequency domain observation signal X (τ, w) send into the structure best submatrix module.

Step 300: the best submatrix module calculating of utilization structure time-frequency domain observation signal vector X (τ, w) the diagonal angle rate and the search largest of correspondence make up best submatrix according to maximum diagonal angle rate; From hybrid matrix A, increase progressively and a mutually different principles of selected M column vector at first successively by the row mark; Construct

individual M * M and tie up submatrix, then each submatrix is carried out respective algorithms and handle; Specifically comprise substep 310,320,330,340,350.

Step 310: at first to first submatrix A of hybrid matrix A ₁Handle; Comprise substep 311,312,313,314.

Step 311: first submatrix A of structure hybrid matrix A ₁, and ask A ₁Inverse matrix

At first, utilize the column vector of hybrid matrix A, mark increases progressively and a mutually different principles of selected M column vector is constructed suc as formula first M shown in (2) * M dimension submatrix A according to being listed as ₁:

A wherein _1i, a _2i..., a _MiBe the i of the hybrid matrix A (column element of 1≤i≤M).Submatrix A ₁Column vector respectively corresponding hybrid matrix A the 1st, 2 ..., the M column element.

Then, ask submatrix A ₁Inverse matrix , because hybrid matrix A is for the row non-singular matrix, so submatrix A ₁There is inverse matrix:

$A_1^{-1}=\mathrm{inv}\left(A_1\right)---\left(3\right)$

Wherein operational character inv () expression finding the inverse matrix operation.

Step 312: utilize step 311 gained matrix

Calculate first and just estimate signal vector S ₁(τ, ω).Just estimate signal vector S ₁(τ ω) is expressed as:

$S_1(\mathrm{\τ},\mathrm{\ω})=A_1^{-1}X(\mathrm{\τ},\mathrm{\ω})---\left(4\right)$

Wherein X (τ, ω)=[X ₁(τ, ω), X ₂(τ, ω) ..., X _M(τ, ω)] ^TThe M that receives for sensor ties up time-frequency domain observation signal vector.Just estimate signal vector S ₁(τ ω) is similarly the M dimensional vector.

Step 313: calculate and just estimate signal vector S ₁(τ, covariance matrix C ω) ₁, shown in (5):

$C_1=S_1{S}_1^H=\left(A_1^{-1}X\right){\left(A_1^{-1}X\right)}^H=A_1^{-1}{\mathrm{XX}}^H{\left(A_1^{-1}\right)}^H---\left(5\right)$

Here, just estimate signal vector S ₁(τ ω) is abbreviated as S ₁(τ w) is abbreviated as X to observation signal vector X; Operational character () ^HThe operation of expression conjugate transpose.

Step 314: calculate covariance matrix C ₁Diagonal angle rate d ₁, shown in (6), (7), (8).

At first, to covariance matrix C ₁Each element get real part, structure covariance matrix C ₁Real part matrix R ₁, shown in (6):

R ₁＝real(C ₁) (6)

Wherein each element real part operation of matrix is got in operational character real () expression.

Then, structure real part matrix R ₁The absolute value matrix M ₁, shown in (7):

M ₁＝abs(R ₁) (7)

Wherein each element absolute value operation of matrix is asked in operational character abs () expression.

At last, by the absolute value matrix M ₁Calculate diagonal angle rate d ₁, shown in (8):

$d_1=\frac{\mathrm{sum}\left(\mathrm{diag}\left(M_1\right)\right)}{\mathrm{sum}\left(M_1\right)}---\left(8\right)$

Wherein all elements and operation in matrix or the vector are asked in operational character sum () expression.The operation of matrix principal diagonal element is got in operational character diag () expression, the column vector of its result for being made up of each main diagonal element of matrix.

Step 320: secondly second sub-matrices of hybrid matrix A is handled; Comprise substep 321,322,323,324.

Step 321: step 311 roughly the same, utilize the column vector of hybrid matrix A, increase progressively and a mutually different principles of selected M column vector is constructed second M * M and tieed up submatrix A according to the row mark ₂, and ask A ₂Inverse matrix

At first, utilize the element among the hybrid matrix A, construct the second sub-matrices A ₂, shown in (9):

A wherein _1i, a _2i..., a _MiI (1≤i≤M+1 and i ≠ M) column element, and A for hybrid matrix A ₂Column vector respectively corresponding hybrid matrix A the 1st, 2 ..., M-1, M+1 column element.

Then, ask submatrix A ₂Inverse matrix because hybrid matrix A for the row non-singular matrix, so submatrix A ₂Be non-singular matrix, have inverse matrix

Shown in (10):

$A_2^{-1}=\mathrm{inv}\left(A_2\right)---\left(10\right)$

Wherein operational character inv () is the finding the inverse matrix operation.

Step 322: utilize matrix

Estimate signal vector S at the beginning of calculating second ₂(τ, ω), shown in (11):

$S_2(\mathrm{\τ},\mathrm{\ω})=A_2^{-1}X(\mathrm{\τ},\mathrm{\ω})---\left(11\right)$

Wherein X (τ, ω)=[X ₁(τ, ω), X ₂(τ, ω) ..., X _M(τ, ω)] ^TThe M that receives for sensor ties up the observation signal vector (with X in the step 312 (τ ω) is same time-frequency domain observation signal vector).Just estimate signal vector S ₂(τ ω) is similarly the M dimensional vector.

Step 323: calculate and just estimate signal vector S ₂(τ, covariance matrix C ω) ₂, shown in (12):

${\mathrm{<figure-callout\; id="2"\; label="C"\; filenames="120223103459.png"\; state="\{\{state\}\}">C</figure-callout>}}_2=S_2{S}_2^H=\left(A_2^{-1}X\right){\left(A_2^{-1}X\right)}^H=A_2^{-1}{\mathrm{XX}}^T{\left(A_2^{-1}\right)}^H---\left(12\right)$

Just estimate signal vector S in the formula ₂(τ, ω) (τ w) simplifies equally, and operational character () with observation signal vector X ^HThe operation of same expression conjugate transpose.

Step 324: calculate covariance matrix C ₂Diagonal angle rate d ₂, shown in (13), (14), (15):

At first, to covariance matrix C ₂Each element get real part, structure covariance matrix C ₂Real part matrix R ₂, shown in (13):

R ₂＝real(C ₂) (13)

Wherein each element real part operation of matrix is got in operational character real () expression.

Secondly, structure real part matrix R ₂The absolute value matrix M ₂, shown in (14):

M ₂＝abs(R ₂) (14)

Wherein each element absolute value operation of matrix is asked in operational character abs () expression.

At last, by the absolute value matrix M ₂Calculate diagonal angle rate d ₂, shown in (15):

${\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="120223103459.png"\; state="\{\{state\}\}">d</figure-callout>}}_2=\frac{\mathrm{sum}\left(\mathrm{diag}\left(M_2\right)\right)}{\mathrm{sum}\left(M_2\right)}---\left(15\right)$

Wherein all elements and operation in matrix or the vector are asked in operational character sum () expression, and the operation of matrix principal diagonal element is got in operational character diag () expression, the column vector of its result for being made up of each main diagonal element of matrix.

Step 330: its complementary submatrix to hybrid matrix A carries out same treatment until

sub-matrices; Processing procedure comprises substep 331,332,333,334.

Step 331: roughly the same step 311, step 321, utilize the column vector of hybrid matrix A, increase progressively and a mutually different principles of selected M column vector constructs the according to the row mark

Individual M * M ties up submatrix, here with p sub-matrices A _pFor example describes, and ask its inverse matrix

, shown in (16), (17).

At first, utilize the element of hybrid matrix A, construct p sub-matrices A _p, shown in (16):

A wherein _1i, a _2i..., a _MiBe the i of the hybrid matrix A (column element of N-M+1≤i≤N)., and submatrix A _pThe column vector N-M+1 of corresponding hybrid matrix A respectively, N-M+2 ..., the N column element.

Then, ask submatrix A _pInverse matrix because hybrid matrix A for the row non-singular matrix, so submatrix A _pBe non-singular matrix, have inverse matrix

Shown in (17):

$A_p^{-1}=\mathrm{inv}\left(A_p\right)---\left(17\right)$

Wherein operational character inv () is the finding the inverse matrix operation.

Step 332: utilize inverse matrix Estimate signal vector S at the beginning of calculating p _p(τ, ω), shown in (18):

$S_p(\mathrm{\&tau;},\mathrm{\&omega;})=A_p^{-1}X(\mathrm{\&tau;},\mathrm{\&omega;})---\left(18\right)$

Wherein X (τ, ω)=[X ₁(τ, ω), X ₂(τ, ω) ..., X _M(τ, ω)] ^TThe M that receives for sensor ties up the observation signal vector (with X in step 312, the step 322 (τ ω) is same time-frequency domain observation signal vector).Just estimate signal vector S _p(τ ω) is similarly the M dimensional vector.

Step 333: calculate and just estimate signal vector S _p(τ, covariance matrix C ω) _p, shown in (19):

$C_p=S_p{S}_p^H=\left(A_p^{-1}X\right){\left(A_p^{-1}X\right)}^H=A_p^{-1}{\mathrm{XX}}^T{\left(A_p^{-1}\right)}^H---\left(19\right)$

Just estimate signal vector S in the formula _p(τ, ω) (τ w) simplifies equally, and operational character () with observation signal vector X ^HThe operation of same expression conjugate transpose.

Step 334: calculate covariance matrix C _pDiagonal angle rate d _p, shown in (20), (21), (22):

At first, to covariance matrix C _pEach element get real part, structure covariance matrix C _pReal part matrix R _p, shown in (20):

R _p＝real(C _p) (20)

Wherein each element real part operation of matrix is got in operational character real () expression.

Secondly, structure real part matrix R _pThe absolute value matrix M _p, shown in (21):

M _p＝abs(R _p) (21)

Wherein each element absolute value operation of matrix is asked in operational character abs () expression.

At last, by the absolute value matrix M _pCalculate diagonal angle rate d _p, shown in (22):

$d_p=\frac{\mathrm{sum}\left(\mathrm{diag}\left(M_p\right)\right)}{\mathrm{sum}\left(M_p\right)}---\left(22\right)$

Wherein all elements and operation in matrix or the vector are asked in operational character sum () expression, and the operation of matrix principal diagonal element is got in operational character diag () expression, the column vector of its result for being made up of each main diagonal element of matrix.

Step 340: 310-330 obtains through step

Individual different diagonal angle rate d ₁, d ₂..., d _p, this step is retrieved maximum diagonal angle rate, and writes down the subscript value of maximum diagonal angle rate, shown in (23):

$\mathrm{index}=\underset{i}{\arg \max }\left\{d_i\right\}(i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,C_N^M)---\left(23\right)$

Wherein index writes down the subscript value of maximum diagonal angle rate.

Step 350: combine maximum diagonal angle rate subscript value index, construct best submatrix A _Opt

Said best submatrix is that subscript value is the submatrix of index, is the index sub-matrices of hybrid matrix A, shown in (24):

A _opt＝A _index (24)

Step 400: utilize the time-frequency domain separation module to carry out the separation of source signal, comprise substep 401,402 at time-frequency domain.

Step 401: utilize best submatrix A _OptMake up time-frequency domain initially-separate signal:

In conjunction with best submatrix A _OptAnd time-frequency domain observation signal vector X (τ, ω), structure time-frequency domain initially-separate signal vector

, shown in (25):

$\tilde{S}(\mathrm{\&tau;},\mathrm{\&omega;}){=A}_{\mathrm{opt}}^{-1}X(\mathrm{\&tau;},\mathrm{\&omega;})---\left(25\right)$

Time-frequency domain initially-separate signal vector wherein $\tilde{S}(\mathrm{\&tau;},\mathrm{\&omega;})={[{\tilde{S}}_1(\mathrm{\&tau;},\mathrm{\&omega;}),{\tilde{S}}_2(\mathrm{\&tau;},\mathrm{\&omega;}),\&CenterDot;\&CenterDot;\&CenterDot;,{\tilde{S}}_M(\mathrm{\&tau;},\mathrm{\&omega;})]}^T$ Be the M dimensional vector, and time-frequency domain source signal vector S (τ, ω)=[S ₁(τ, ω), S ₂(τ, ω) ..., S _N(τ, ω)] ^TBe the N dimensional vector.So there is the assignment problem between time-frequency domain initially-separate signal vector and the final separation signal vector of time-frequency domain.

Step 402: utilize the time-frequency domain distribution module, solve the assignment problem of the final separation signal vector of time-frequency domain initially-separate signal vector and time-frequency domain.Distribution principle is following:

With best submatrix A _OptThe row that correspond among the hybrid matrix A of column vector be labeled as index1, index2 ..., indexM, right back-pushed-type (26) distributes:

${\hat{S}}_{\mathrm{indexi}}(\mathrm{\&tau;},\mathrm{\&omega;})={\tilde{S}}_i(\mathrm{\&tau;},\mathrm{\&omega;})---\left(26\right)$

I=1 wherein, 2 ..., M,

The final separation signal vector of expression time-frequency domain $\hat{S}(\mathrm{\&tau;},\mathrm{\&omega;})=[{\hat{S}}_1(\mathrm{\&tau;},\mathrm{\&omega;}),{\hat{S}}_2(\mathrm{\&tau;},\mathrm{\&omega;}),\&CenterDot;\&CenterDot;\&CenterDot;,{\hat{S}}_N(\mathrm{\&tau;},\mathrm{\&omega;})]$ Indexi element,

Expression

I element.Formula (26) is intended to indexi element with

and is changed to i the element of

and equates that all the other elements put 0.

Step 500: utilize the source signal estimation module to obtain the estimation of each road source signal, and export.

Circulation execution in step 300, step 400 until all the time frequency dispose; The final separation signal Vector Groups of all time-frequency domains is combined into the time-frequency domain separation matrix; Its row vector is the time-frequency domain separation signal; Be designated as

the time-frequency domain separation signal passed through Fourier inversion in short-term, obtain being estimated as of each road source signal:

${\hat{s}}_j=\mathrm{ISTFT}\left[{\hat{S}}_j(\mathrm{\&tau;},\mathrm{\&omega;})\right],j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N---\left(27\right)$

Wherein

is j separating obtained time-frequency domain signal, and

is the time domain separation signal.

Beneficial effect: frequency was the stack of which road source signal when the present invention can know each exactly, had reduced the requirement to source signal statistical property and sparse property, was correlated with at source signal and perhaps also can separates preferably under the not strict sparse situation.The present invention is equally applicable to the very little situation of angle between the hybrid matrix column vector.In addition, the processing of the present invention frequency and each sub-matrices during to each can executed in parallel, is beneficial to hardware and realizes.

Description of drawings

Fig. 1 is the blind source separation method synoptic diagram of instantaneous mixing.

Fig. 2 is the present invention owes blind separation surely to instantaneous mixed signal a system chart.

Fig. 3 is the blind surely separation emulation experiment figure that owes of four road voice signal linear instantaneous mixing.

Wherein, Fig. 3 (a) four road source signals; Fig. 3 (b) three road observation signals; Fig. 3 (c) four tunnel separation signals.

Fig. 4 is the blind surely separation emulation experiment figure that owes of four road correlated source signal linear instantaneous mixing.

Wherein, Fig. 4 (a) four road source signals; Fig. 4 (b) three road observation signals; Fig. 4 (c) four tunnel separation signals.

Fig. 5 is the blind surely emulation experiment figure that separates that owes of one road white Gaussian noise signal and the mixing of three road voice signal linear instantaneous.

Wherein, Fig. 5 (a) four road source signals; Fig. 5 (b) three road observation signals; Fig. 5 (c) four tunnel separation signals.

Fig. 6 is the blind surely emulation experiment figure that separates that owes of two-way music signal and the mixing of two-way voice signal linear instantaneous.

Wherein, Fig. 6 (a) four road source signals; Fig. 6 (b) three road observation signals; Fig. 6 (c) four tunnel separation signals.

Fig. 7 is the emulation experiment figure to the situation that angle is all very little between all column vectors of hybrid matrix.

Wherein, Fig. 7 (a) three road source signals; Fig. 7 (b) two road observation signals; Fig. 7 (c) three tunnel separation signals.

Embodiment

Below in conjunction with accompanying drawing and embodiment the present invention is elaborated.

Basic thought of the present invention:, observation signal is carried out Short Time Fourier Transform be equivalent to and earlier source signal carried out Short Time Fourier Transform premultiplication hybrid matrix again according to the linear behavio(u)r of Short Time Fourier Transform.Based on this linear behavio(u)r, each the time frequency can regard the linear superposition (if some time frequency only contain one road source signal, then can regard the weighting of this source signal as) of several roads source signal as.Frequency was the stack (comprising the situation that only contains one road source signal) of which road source signal but also the time-frequency value that can estimate source signal when the present invention not only can judge each exactly.

The present invention at first transforms to time-frequency domain with observation signal; Secondly increase progressively according to the row mark and individual M * M that a mutually different principles of selected M column vector is constructed hybrid matrix ties up submatrix, and ask the inverse matrix of each sub-matrices.Estimate signal vector at the beginning of then each inverse matrix and observation signal multiplication of vectors to be constructed

individual and ask each just to estimate covariance matrix and real part matrix and absolute value matrix of signal vector.The covariance matrix of the present invention's hypothesis equals just to estimate the product of signal vector and its conjugate transpose.The real part matrix is the matrix that each element real part of covariance matrix constitutes, and the absolute value matrix is the matrix that each element absolute value of real part matrix constitutes.Obtain individual diagonal angle rate by individual absolute value matrix computations at last; Get the corresponding submatrix of maximum diagonal angle rate and just estimate signal vector, and M is tieed up the initially-separate signal vector be assigned as N to tie up final separation signal vectorial respectively as best submatrix and initially-separate signal vector.

The present invention proposes, and to owe the system chart of blind separation surely to instantaneous mixed signal as shown in Figure 2; This method basic procedure is following: from sound bank, choose different audio signals as source signal; Set up hybrid matrix A on their own, in time domain source signal is carried out the linear instantaneous mixing and obtain observation signal.The time domain observation signal is arrived time-frequency domain through Short Time Fourier Transform, and mark increases progressively and a mutually different principles of selected M column vector is constructed the submatrix and the inverse matrix thereof of hybrid matrix according to being listed as afterwards, and asks and just estimate signal vector.Ask each covariance matrix of just estimating signal vector and real part matrix, absolute value matrix then, and then ask diagonal angle rate, best submatrix, initially-separate signal vector, and the initially-separate signal vector is assigned as final separation signal vector.At last the time-frequency domain separation signal is obtained the estimation of source signal through Fourier inversion in short-term to time domain.

In order to verify that the present invention proposes the validity of method, has carried out computer simulation experiment.Emulation experiment is mixed into three the tunnel and three the tunnel with four road signals, and to be mixed into two the tunnel be example, and source signal not only is included in the sparse property of time-frequency domain voice signal preferably, and also have noise and music signal, even comprise the correlated source signal.

Embodiment 1

The voice signal of from sound bank, choosing 4 tunnel different speakers is got 50000 points respectively as source signal, carries out linear instantaneous in time domain and mixes, and hybrid matrix is following:

$A=\left[\begin{array}{cccc}0.7875& 0.6368& 0.6117& 0.4277\\ 0.4161& 0.4792& 0.7146& 0.7530\\ 0.4545& 0.6039& 0.3391& 0.5000\end{array}\right]$

The absolute value of the determinant of each sub-matrices of hybrid matrix A is respectively 0.0744,0.0697,0.0838,0.0815.Three tunnel mixed signals that obtain are arrived time-frequency domain through Short Time Fourier Transform, and the time-domain signal frame length of Short Time Fourier Transform is got 1024 points, overlapping 768 points of interframe, and frequency domain resolution is got 1024 points.Four road original voice signal time domain waveforms are shown in accompanying drawing 3 (a), and the time domain waveform of carrying out mixed three road observation signals of linear instantaneous is shown in accompanying drawing 3 (b), and isolated four road voice signal time domain waveforms of the present invention are shown in accompanying drawing 3 (c).Comparison diagram 3 (a) and Fig. 3 (c) can find out that restorability is satisfactory.In addition, the present invention also utilizes output signal-to-noise ratio that the separating property of this method has been carried out objective measurement, and the computing formula of output signal-to-noise ratio is suc as formula shown in (21):

${\mathrm{SNR}}_i=10{\log }_{10}\left(\frac{\underset{n=1}{\overset{N_S}{\mathrm{\&Sigma;}}}{\left|s_i\left(n\right)\right|}^2}{\underset{n=1}{\overset{N_S}{\mathrm{\&Sigma;}}}{|s_i\left(n\right)-{\tilde{s}}_i\left(n\right)|}^2}\right),i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N---\left(29\right)$

In the formula, s _i(n) be i original signal,

The estimation of i the original signal that obtains for blind separation.N _SFor signals sampling is counted, N is the source signal number.For fear of the uncertain influence that brings of amplitude, s here _i(n) and

Carrying out all having passed through average power normalization processing before SNR calculates.The output signal-to-noise ratio of separation signal is as shown in table 1.

Table 1 four tunnel different phonetic output signal-to-noise ratio SNR

Embodiment 2

One of them advantage of the present invention is to have reduced the requirement to the source signal statistical property, can owe to decide separation to the correlated source signal.From sound bank, choose one road voice signal, the different time sections of getting these voice then constitutes four road correlated source signals.Four road original coherent signal time domain waveforms are shown in accompanying drawing 4 (a), and the time domain waveform of carrying out mixed three road observation signals of linear instantaneous is shown in accompanying drawing 4 (b), and isolated four road time domain plethysmographic signals of the present invention are shown in accompanying drawing 4 (c).The output signal-to-noise ratio of separation signal is as shown in table 2.

Table 2 four tunnel related voice output signal-to-noise ratio SNR

Embodiment 3

What the present invention can solve weak sparse signal preferably owes to decide separation.Be that example is explained advantage of the present invention with the relatively poor white Gaussian noise of sparse property below.From sound bank, choose one road white Gaussian noise and three tunnel different phonetic signals.A road original white Gaussian noise and three road voice signal time domain waveforms are shown in accompanying drawing 5 (a); The time domain waveform of carrying out mixed three road observation signals of linear instantaneous is shown in accompanying drawing 5 (b), and isolated four road time domain plethysmographic signal figure of the present invention are shown in accompanying drawing 5 (c).The output signal-to-noise ratio of separation signal is as shown in table 3.

Table 3 one road noise, three tunnel voice output signal to noise ratio snr

Embodiment 4

Be that example further specifies the advantage that the present invention is directed to weak sparse signal with the relatively poor music signal of sparse property below.From sound bank, choose different music signals of two-way and two-way different phonetic signal.Original two-way music signal and two-way voice signal time domain waveform are shown in accompanying drawing 6 (a), and the time domain waveform of carrying out mixed three road observation signals of linear instantaneous is shown in accompanying drawing 6 (b), and isolated four road time domain plethysmographic signals of the present invention are shown in accompanying drawing 6 (c).The output signal-to-noise ratio of separation signal is as shown in table 4.

Table 4 two-way music, two-way voice output signal to noise ratio snr

Embodiment 5

Be the validity that example explanation the present invention is directed to this type of situation with the hybrid matrix that angle is all very little between all column vectors below.Hybrid matrix is following:

$A=\left[\begin{array}{ccc}\cos (30\&times;\mathrm{\&pi;}/180)& \cos (31\&times;\mathrm{\&pi;}/180)& \cos (32\&times;\mathrm{\&pi;}/180)\\ \sin (30\&times;\mathrm{\&pi;}/180)& \sin (31\&times;\mathrm{\&pi;}/180)& \sin (32\&times;\mathrm{\&pi;}/180)\end{array}\right]$

Among the hybrid matrix A between first column vector and secondary series vector angle be 1 degree, angle is similarly 1 and spends between secondary series vector and the 3rd column vector.The absolute value of the determinant of each sub-matrices is respectively 0.0175,0.0349,0.0175.From sound bank, choose three tunnel different phonetic signals.Original voice signal time domain waveform is shown in accompanying drawing 7 (a), and the time domain waveform of carrying out mixed two road observation signals of linear instantaneous is shown in accompanying drawing 7 (b), and isolated three road time domain plethysmographic signals of the present invention are shown in accompanying drawing 7 (c).The output signal-to-noise ratio of separation signal is as shown in table 5.

Table 5 clip angle hybrid matrix output signal-to-noise ratio SNR

Shown by above analysis and experimental data: the present invention is applicable to that linear instantaneous mixes and owes blind surely separation; The situation that can under signal correction or not strict sparse situation, for example add one road noise or two-way music is separated preferably, and can handle all very little situation of angle between all column vectors of hybrid matrix.

Sum up: the present invention is that the instantaneous mixing with time domain is transformed into time-frequency domain and handles, and on the basis that hypothesis hybrid matrix A has estimated to obtain through classic method, the present invention is through the submatrix and the inverse matrix thereof of structure hybrid matrix, and asks and just estimate signal vector.Ask each covariance matrix of just estimating signal vector and real part matrix, absolute value matrix then, and then ask maximum diagonal angle rate.Get the corresponding submatrix of maximum diagonal angle rate and just estimate signal vector, and M is tieed up the initially-separate signal vector be assigned as N to tie up final separation signal vectorial respectively as best submatrix and initially-separate signal vector.The present invention has reduced the requirement to source signal statistical property and sparse property, is correlated with at source signal and perhaps also can separates preferably under the not strict sparse situation.The present invention is equally applicable to the very little situation of angle between the hybrid matrix column vector.In addition, the processing of the present invention frequency and each sub-matrices during to each can executed in parallel, is beneficial to hardware and realizes.

Above content is to combine optimal technical scheme to the further explain that the present invention did, and can not assert that the practical implementation of invention only limits to these explanations.Under the present invention, the those of ordinary skill of technical field, under the prerequisite that does not break away from design of the present invention, simple deduction and replacement can also be made, all protection scope of the present invention should be regarded as.

Claims (1)

1. One kind based on maximum matrix diagonal angle rate owe to decide blind separating method, may further comprise the steps:

   Step 100: utilize the observation signal of the instantaneous mixing of sensor, be expressed as:

   x(t)＝As(t) (1)

   X (t) in the formula (1)=[x ₁(t), x ₂(t) ..., x _M(t)] ^TThe M dimension observation signal vector that expression is received by sensor, A are represented M * N (the normalized capable full rank hybrid matrix of dimension of M＜N), s (t)=[s ₁(t), s ₂(t) ..., s _N(t)] ^TThe source signal vector that expression N dimension is unknown;

   Step 200; The observation signal of instantaneous mixing is sent into the orthogonal transformation module carries out Short Time Fourier Transform, obtain time-frequency domain observation signal X (τ, w); And with this time-frequency domain observation signal X (τ, w) send into the structure best submatrix module;

   Step 300: the best submatrix module calculating of utilization structure time-frequency domain observation signal vector X (τ, w) the diagonal angle rate and the search largest of correspondence make up best submatrix according to maximum diagonal angle rate; From hybrid matrix A, increase progressively and a mutually different principles of selected M column vector at first successively by the row mark; Construct

   

   individual M * M and tie up submatrix, then each submatrix is carried out respective algorithms and handle; Specifically comprise substep 310,320,330,340,350;

   Step 310: at first to first submatrix A of hybrid matrix A ₁Handle; Comprise substep 311,312,313,314:

   Step 311: first submatrix A of structure hybrid matrix A ₁, and ask A ₁Inverse matrix

   At first, utilize the column vector of hybrid matrix A, mark increases progressively and a mutually different principles of selected M column vector is constructed suc as formula first M shown in (2) * M dimension submatrix A according to being listed as ₁:

   

   A wherein _1i, a _2i..., a _MiBe the i of the hybrid matrix A (column element of 1≤i≤M); Submatrix A ₁Column vector respectively corresponding hybrid matrix A the 1st, 2 ..., the M column element;

   Then, ask submatrix A ₁Inverse matrix

   

   Because hybrid matrix A is for the row non-singular matrix, so submatrix A ₁There is inverse matrix:

   $A_1^{-1}=\mathrm{inv}\left(A_1\right)---\left(3\right)$

   Wherein operational character inv () expression finding the inverse matrix operation;

   Step 312: utilize step 311 gained matrix Calculate first and just estimate signal vector S ₁(τ, ω); Just estimate signal vector S ₁(τ ω) is expressed as:

   $S_1(\mathrm{\&tau;},\mathrm{\&omega;})=A_1^{-1}X(\mathrm{\&tau;},\mathrm{\&omega;})---\left(4\right)$

   Wherein X (τ, ω)=[X ₁(τ, ω), X ₂(τ, ω) ..., X _M(τ, ω)] ^TThe M that receives for sensor ties up time-frequency domain observation signal vector; Just estimate signal vector S ₁(τ ω) is similarly the M dimensional vector;

   Step 313: calculate and just estimate signal vector S ₁(τ, covariance matrix C ω) ₁, shown in (5):

   $C_1=S_1{S}_1^H=\left(A_1^{-1}X\right){\left(A_1^{-1}X\right)}^H=A_1^{-1}{\mathrm{XX}}^H{\left(A_1^{-1}\right)}^H---\left(5\right)$

   Here, just estimate signal vector S ₁(τ ω) is abbreviated as S ₁(τ w) is abbreviated as X to observation signal vector X; Operational character () ^HThe operation of expression conjugate transpose;

   Step 314: calculate covariance matrix C ₁Diagonal angle rate d ₁, shown in (6), (7), (8):

   At first, to covariance matrix C ₁Each element get real part, structure covariance matrix C ₁Real part matrix R ₁, shown in (6):

   R ₁＝real(C ₁) (6)

   Wherein each element real part operation of matrix is got in operational character real () expression;

   Then, structure real part matrix R ₁The absolute value matrix M ₁, shown in (7):

   M ₁＝abs(R ₁) (7)

   Wherein each element absolute value operation of matrix is asked in operational character abs () expression;

   At last, by the absolute value matrix M ₁Calculate diagonal angle rate d ₁, shown in (8):

   $d_1=\frac{\mathrm{sum}\left(\mathrm{diag}\left(M_1\right)\right)}{\mathrm{sum}\left(M_1\right)}---\left(8\right)$

   Wherein all elements and operation in matrix or the vector are asked in operational character sum () expression; The operation of matrix principal diagonal element is got in operational character diag () expression, the column vector of its result for being made up of each main diagonal element of matrix;

   Step 320: secondly second sub-matrices of hybrid matrix A is handled; Comprise substep 321,322,323,324:

   Step 321: step 311 roughly the same, utilize the column vector of hybrid matrix A, increase progressively and a mutually different principles of selected M column vector is constructed second M * M and tieed up submatrix A according to the row mark ₂, and ask A ₂Inverse matrix

   

   At first, utilize the element among the hybrid matrix A, construct the second sub-matrices A ₂, shown in (9):

   

   A wherein _1i, a _2i..., a _MiI (1≤i≤M+1 and i ≠ M) column element, and A for hybrid matrix A ₂Column vector respectively corresponding hybrid matrix A the 1st, 2 ..., M-1, M+1 column element;

   Then, ask submatrix A ₂Inverse matrix because hybrid matrix A for the row non-singular matrix, so submatrix A ₂Be non-singular matrix, have inverse matrix Shown in (10):

   $A_2^{-1}=\mathrm{inv}\left(A_2\right)---\left(10\right)$

   Wherein operational character inv () is the finding the inverse matrix operation;

   Step 322: utilize matrix

   

   Estimate signal vector S at the beginning of calculating second ₂(τ, ω), shown in (11):

   $S_2(\mathrm{\&tau;},\mathrm{\&omega;})=A_2^{-1}X(\mathrm{\&tau;},\mathrm{\&omega;})---\left(11\right)$

   Wherein X (τ, ω)=[X ₁(τ, ω), X ₂(τ, ω) ..., X _M(τ, ω)] ^TThe M that receives for sensor ties up the observation signal vector (with X in the step 312 (τ ω) is same time-frequency domain observation signal vector); Just estimate signal vector S ₂(τ ω) is similarly the M dimensional vector;

   Step 323: calculate and just estimate signal vector S ₂(τ, covariance matrix C ω) ₂, shown in (12):

   $C_2=S_2{S}_2^H=\left(A_2^{-1}X\right){\left(A_2^{-1}X\right)}^H=A_2^{-1}{\mathrm{XX}}^T{\left(A_2^{-1}\right)}^H---\left(12\right)$

   Just estimate signal vector S in the formula ₂(τ, ω) (τ w) simplifies equally, and operational character () with observation signal vector X ^HThe operation of same expression conjugate transpose;

   Step 324: calculate covariance matrix C ₂Diagonal angle rate d ₂, shown in (13), (14), (15):

   At first, to covariance matrix C ₂Each element get real part, structure covariance matrix C ₂Real part matrix R ₂, shown in (13):

   R ₂＝real(C ₂) (13)

   Wherein each element real part operation of matrix is got in operational character real () expression;

   Secondly, structure real part matrix R ₂The absolute value matrix M ₂, shown in (14):

   M ₂＝abs(R ₂) (14)

   Wherein each element absolute value operation of matrix is asked in operational character abs () expression;

   At last, by the absolute value matrix M ₂Calculate diagonal angle rate d ₂, shown in (15):

   $d_2=\frac{\mathrm{sum}\left(\mathrm{diag}\left(M_2\right)\right)}{\mathrm{sum}\left(M_2\right)}---\left(15\right)$

   Wherein all elements and operation in matrix or the vector are asked in operational character sum () expression, and the operation of matrix principal diagonal element is got in operational character diag () expression, the column vector of its result for being made up of each main diagonal element of matrix;

   Step 330: its complementary submatrix to hybrid matrix A carries out same treatment until

   

   sub-matrices; Processing procedure comprises substep 331,332,333,334:

   Step 331: roughly the same step 311, step 321, utilize the column vector of hybrid matrix A, increase progressively and a mutually different principles of selected M column vector constructs the according to the row mark

   

   Individual M * M ties up submatrix, here with p sub-matrices A _pFor example describes, and ask its inverse matrix Shown in (16), (17):

   At first, utilize the element of hybrid matrix A, construct p sub-matrices A _p, shown in (16):

   

   A wherein _1i, a _2i..., a _MiI (N-M+1≤i≤N) column element, and submatrix A for hybrid matrix A _pThe column vector N-M+1 of corresponding hybrid matrix A respectively, N-M+2 ..., the N column element;

   Then, ask submatrix A _pInverse matrix because hybrid matrix A for the row non-singular matrix, so submatrix A _pBe non-singular matrix, have inverse matrix

   

   Shown in (17):

   $A_p^{-1}=\mathrm{inv}\left(A_p\right)---\left(17\right)$

   Wherein operational character inv () is the finding the inverse matrix operation;

   Step 332: utilize inverse matrix

   

   Estimate signal vector S at the beginning of calculating p _p(τ, ω), shown in (18):

   $S_p(\mathrm{\&tau;},\mathrm{\&omega;})=A_p^{-1}X(\mathrm{\&tau;},\mathrm{\&omega;})---\left(18\right)$

   Wherein X (τ, ω)=[X ₁(τ, ω), X ₂(τ, ω) ..., X _M(τ, ω)] ^TThe M that receives for sensor ties up the observation signal vector (with X in step 312, the step 322 (τ ω) is same time-frequency domain observation signal vector); Just estimate signal vector S _p(τ ω) is similarly the M dimensional vector;

   Step 333: calculate and just estimate signal vector S _p(τ, covariance matrix C ω) _p, shown in (19):

   $C_p=S_p{S}_p^H=\left(A_p^{-1}X\right){\left(A_p^{-1}X\right)}^H=A_p^{-1}{\mathrm{XX}}^T{\left(A_p^{-1}\right)}^H---\left(19\right)$

   Just estimate signal vector S in the formula _p(τ, ω) (τ w) simplifies equally, and operational character () with observation signal vector X ^HThe operation of same expression conjugate transpose;

   Step 334: calculate covariance matrix C _pDiagonal angle rate d _p, shown in (20), (21), (22):

   At first, to covariance matrix C _pEach element get real part, structure covariance matrix C _pReal part matrix R _p, shown in (20):

   R _p＝real(C _p) (20)

   Wherein each element real part operation of matrix is got in operational character real () expression;

   Secondly, structure real part matrix R _pThe absolute value matrix M _p, shown in (21):

   M _p＝abs(R _p) (21)

   Wherein each element absolute value operation of matrix is asked in operational character abs () expression;

   At last, by the absolute value matrix M _pCalculate diagonal angle rate d _p, shown in (22):

   $d_p=\frac{\mathrm{sum}\left(\mathrm{diag}\left(M_p\right)\right)}{\mathrm{sum}\left(M_p\right)}---\left(22\right)$

   Wherein all elements and operation in matrix or the vector are asked in operational character sum () expression, and the operation of matrix principal diagonal element is got in operational character diag () expression, the column vector of its result for being made up of each main diagonal element of matrix;

   Step 340: 310-330 obtains through step

   

   Individual different diagonal angle rate d ₁, d ₂..., d _p, this step is retrieved maximum diagonal angle rate, and writes down the subscript value of maximum diagonal angle rate, shown in (23):

   $\mathrm{index}=\underset{i}{\arg \max }\left\{d_i\right\}(i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,C_N^M)---\left(23\right)$

   Wherein index writes down the subscript value of maximum diagonal angle rate;

   Step 350: combine maximum diagonal angle rate subscript value index, construct best submatrix A _Opt

   Said best submatrix is that subscript value is the submatrix of index, is the index sub-matrices of hybrid matrix A, shown in (24):

   A _opt＝A _index (24)

   Step 400: utilize the time-frequency domain separation module to carry out the separation of source signal, comprise substep 401,402 at time-frequency domain;

   Step 401: utilize best submatrix A _OptMake up time-frequency domain initially-separate signal:

   In conjunction with best submatrix A _OptAnd time-frequency domain observation signal vector X (τ, ω), structure time-frequency domain initially-separate signal vector

   

   Shown in (25):

   $\tilde{S}(\mathrm{\&tau;},\mathrm{\&omega;})=A_{\mathrm{opt}}^{-1}X(\mathrm{\&tau;},\mathrm{\&omega;})---\left(25\right)$

   Time-frequency domain initially-separate signal vector wherein $\tilde{S}(\mathrm{\&tau;},\mathrm{\&omega;})={[{\tilde{S}}_1(\mathrm{\&tau;},\mathrm{\&omega;}),{\tilde{S}}_2(\mathrm{\&tau;},\mathrm{\&omega;}),\&CenterDot;\&CenterDot;\&CenterDot;,{\tilde{S}}_M(\mathrm{\&tau;},\mathrm{\&omega;})]}^T$ Be the M dimensional vector, and time-frequency domain source signal vector S (τ, ω)=[S ₁(τ, ω), S ₂(τ, ω) ..., S _N(τ, ω)] ^TBe the N dimensional vector; So there is the assignment problem between time-frequency domain initially-separate signal vector and the final separation signal vector of time-frequency domain;

   Step 402: utilize the time-frequency domain distribution module, solve the assignment problem of the final separation signal vector of time-frequency domain initially-separate signal vector and time-frequency domain; Distribution principle is following:

   With best submatrix A _OptThe row that correspond among the hybrid matrix A of column vector be labeled as index1, index2 ..., indexM, right back-pushed-type (26) distributes:

   ${\hat{S}}_{\mathrm{indexi}}(\mathrm{\&tau;},\mathrm{\&omega;})={\tilde{S}}_i(\mathrm{\&tau;},\mathrm{\&omega;})---\left(26\right)$

   I=1 wherein, 2 ..., M,

   

   The final separation signal vector of expression time-frequency domain $\hat{S}(\mathrm{\&tau;},\mathrm{\&omega;})=[{\hat{S}}_1(\mathrm{\&tau;},\mathrm{\&omega;}),{\hat{S}}_2(\mathrm{\&tau;},\mathrm{\&omega;}),\&CenterDot;\&CenterDot;\&CenterDot;,{\hat{S}}_N(\mathrm{\&tau;},\mathrm{\&omega;})]$ Indexi element,

   

   Expression

   

   I element; Formula (26) be intended to

   

   Indexi element be changed to

   

   I element equate that all the other elements put 0;

   Step 500: utilize the source signal estimation module to obtain the estimation of each road source signal, and export:

   Circulation execution in step 300, step 400 until all the time frequency dispose; The final separation signal Vector Groups of all time-frequency domains is combined into the time-frequency domain separation matrix; Its row vector is the time-frequency domain separation signal; Be designated as

   

   the time-frequency domain separation signal passed through Fourier inversion in short-term, obtain being estimated as of each road source signal:

   ${\hat{s}}_j=\mathrm{ISTFT}\left[{\hat{S}}_j(\mathrm{\&tau;},\mathrm{\&omega;})\right],j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N---\left(27\right)$

   Wherein

   

   is j separating obtained time-frequency domain signal, and

   

   is the time domain separation signal.

Images