Algorithm for source recovery in underdetermined blind source separation based on plane pursuit

FU Weihong,WEI Juan,LIU Naian,and CHEN Jiehu

1.State Key Lab of Integrated Service Networks,Xidian University,Xi’an 710071,China;2.No.10 Research Institution,China Electronics Technology Group Corporation,Chengdu 610036,China

1.Introduction

Blind source separation(BSS)aims to recover sources from a multi-channel system using sensor signals with very little prior knowledge.It has received wide attention for its potential applications in array signal processing[1],speech signal processing[2],biomedical engineering[3],and so on.If the number of sources is larger than that of the sensors,the BSS problem becomes more difficult to solve,and is referred to as the underdetermined blind source separation(UBSS).

Generally,a two-step method,i.e.,estimating the mixing matrix and recovering the source signals,is suggested to solve the UBSS problem.For the first step,many efficient methods have been proposed and achieve good performance[4,5].However,the second step has not yet been well addressed,and the existed algorithms cannot simultaneously achieve high recovery accuracy with low computational complexity.In this paper,we assume that the mixing matrix has been well estimated without estimation error,and focus on the second step,i.e.,recovery of source signals.More specifically,we consider the scenario where the source signal is sparse,i.e.,majorities of the entries in the source signal vector(or its transformed version in a certain domain)are zeros.This scenario is common in reality.In this sparse scenario,UBSS can be commonly regarded as a sparse signal reconstruction problem due to the similar mathematical models between UBSS and compressed sensing(CS)[6].

Currently,there have been many sparse-signal reconstruction algorithms for the UBSS problem.In the early literature,l0norm is proposed to reconstruct the sparse signals,but it turns out to be a nondeterministic polynomial-time hard(NP-hard)problem[7].Since the l0norm can be converted to convex l1norm,the basis pursuit(BP)[8]algorithm is proposed.Although the BP algorithm can obtain the global optimal solution,its computational complexity is huge to meet the requirement of instantaneity.Compared with the BP algorithms,the greedy algorithms[9]have low computational complexity.A typical greedy algorithm is the orthogonal matching pursuit(OMP)algorithm[10].Reconstruction accuracy of OMP is low when the source signals are not sparse enough.Afterwards,more sophisticated greedy algorithms are been developed,such as compressed sampling matching pursuit(CoSaMP)[11]and subspace pursuit(SP)[12].These algorithms have better accuracy performance than the OMP algorithm.However,in these algorithms,the sparsity K should be known in advance,which may not be feasible in many practical applications.Recently,adaptive greedy algorithms are proposed,such as the adaptive sparsity matching pursuit(ASMP)[13]algorithm.Compared with the OMP algorithm,ASMP can guarantee the feasibility of chosen atoms while introducing more computational time.In order to obtain the low complexity and high recovery accuracy,plane pursuit(PP)is proposed in this paper.The PP algorithm selects the atoms according to the correlation between received signals and hyper plane composed by columns of the mixing matrix,which is suitable for the case that the source is not sparse enough and reduces the number of iterations.

2.System model of UBSS

The general objective of BSS is to estimate the source signals from their instantaneous,convolutive,or non-linear measurements which are unknown mixtures of the source signals.Consider the following linear instantaneous mixing system:

where the source signals=[s1(t),s2(t),...,sN(t)]Tis the N-dimentional vector and is assumed to be sparse in time domain or a certain transformation domain.Ais an M×N row-full-rank mixing matrix,x(t)=[x1(t),x2(t),...,xM(t)]Tis the M-dimentionalvector of mixed signals observed by M sensors,and the superscript(.)Tdenotes the transpose operator.

The BSS is called UBSS if the number of sensors M is smaller than that of sources N.The source signal recovery of UBSS aims to recover the source signalssfrom the observedsignalsxand mixing matrixA.Since the mixing matrixAis not of full columnrankin the underdetermined case,the source signalss(t)cannot be obtained through multiplying the mixturesx(t)by the pseudo-inversion ofA.Therefore,recovering the source signals is a big challenge even if the mixing matrixAis known.

3.PP algorithm

Greedy algorithms are widely used in source recovery in UBSS due to low complexity.The typical greedy algorithms,such as the OMP algorithm and the SP algorithm,select matching atoms by calculating the absolute value of the inner product between atoms and the residual.When the source signal is not sparse enough,i.e.,there is more than one nonzero element at each sampling time,we may not select the correct atoms according to the absolute value of the inner product.In order to solve this problem,the ASMP algorithm takes advantage of backtracking to refine the chosen supports and the current approximation in the process.The backtracking process could select more accurate atoms.Unfortunately,this process will dramatically increase the number of iterations. Compared with the OMP algorithm,ASMP has higher reconstruction accuracy but high computational complexity.If we select the atoms accordingto the correlation between sensor signals and hyper plane composed by the column vectors of the mixing matrix,we can improve the convergences peed and guarantee the reliability of the support.Motived by the above,the PP is proposed.

(i)Noiseless case when sources are no more than M-1

Consideringan M×N mixing matrixA,if anyM×M sub-matrix ofAis of full rank,the sub-matrixAi(i=1,2,...,),which is composed by M-1 column vectors ofA,spans a subspace hyper plane Si.ForAi,there exists a unique unit normal vectornisatisfying the following equation:

Assuming that the amount of the nonzero elements of source signals is no more than M-1 and there is no noise at time t,the sensor signalx(t)lies in one or a few of the hyper planes.These hyper planes contain all column vectors,andx(t)is the linear combination of these column vectors.Because the normal vectors of these hyper planes are orthogonal tox(t),we can find out the normal vectors by calculating inner products between thenormal vectors andx(t).If one normal vectornisatis fies

x(t)lies in the ith hyper plane Siand can be expressed as the linear combination of column vectors ofAi:

whereAi= [ai1,ai2,...aiM-1]is the sub-matrix corresponding to normalvectorni.The elements si1(t),si2(t),...,siM-1(t)can be obtained by the least square solution:

The superscript H denotes the conjugate transpose operator.Because the amount of nonzero elements of source signal Ntis no more than M-1,the M-1 elements obtained from(5)include all Ntnonzero elements and M-1-Ntzero elements ins(t).Other N-M+1 elements ofs(t)mustbezeros.DenoteΛ ={i1,i2,...,iM-1}as the subscript set.The M-1 elements of the estimated source signal?s(t)corresponding to Λ are si1(t),si2(t),...,siM-1(t)and the other N-M+1 elements are zeros.

Through above analysis,we can conclude that when the amount of the nonzero elements of source signals is no more than M-1 in the absence of noise,we can accurately identify the hyper plane Siwhich contains all the correct atoms.

(ii)Noisy case when sources are no more than M-1

Now,let us consider the case when the source signal contains no more than M-1 nonzero elements under the condition of noise,and we assume the noise is of Gaussian type.In this case,x(t)would not belong to any hyper plane and we cannot find the normal vectornisatisfying(3).Therefore,we should approximate the source signal by estimating wherex(t)lies in.Assume thatx(t)lies in the ith hyperplane Si.Following(4),x(t)can be expressed as the linear combination of column vectors ofAi:

wherevis the white Gaussian noise with zero mean and variance σ2,and we havex(t)=v.Let li=.Since liis a linear combination of M Gaussian variables,liobeys Gaussian distribution.We can compute E(li)=E(v)=0 and E((-E(li))2)=.We thus can get the probability density function of li:

Because of γ1,(8)can be simplified as

Then the γ1estimated values of source signalscan be calculated by

(iii)The case when sources are more than M-1

In case that the number of nonzero elements is more than γ1,we cannot find the normal vector γ2satisfying(3)no matter the noise exists or not.We should first search for the index i(1)of most matching hyper plane Si(1)by(9).Moreover,the M-1 estimated values of source signals,denoted bycan be obtained by(10).We then calculate the source signalwhich includes M-1 elementsand N-M+1 zero elements.The residual isR1(t)=Because the number of nonzero elements is more than M-1,the residual may not fall below the desired error bound η.We thus identify the index i(2)of the hyper plane,which matches the hyper plane with the residualR1(t)with the highest probability,and we have

After gettingR2(t)=R1(t)-A(2)(t),we then check whether ‖R2‖2＜ η,where η is the desired error bound and is set to 0.01 in the following simulations.The calculation process continues until the norm of the residual falls below η.We finally get the recovery signalwhere P is the total number of iterations.

4.Convergence of the PP algorithm

We now analyze the convergence of the proposed algorithm.Assume that the residual after the kth iteration is represented byRk(t),and the(k+1)th residual is denoted byRk+1(t).For simplicity,we remove t.Then we get

Then

5.Comparison of computational complexity

In this section,we will compare the PP algorithm with OMP,CMP and ASMP algorithms.Considering the computational costs of these four algorithms mainly depend on the floating point multiplication and floating point additions,through analysis we find that the number of floating point multiplications in four algorithms nearly equals the number of additions,we can use the ratio of the number of multiplication to represent the approximate ratio of computational complexity.The computational cost of the PP algorithm primarily depends on Step 2 and Step 3.The number of the multiplication at each iteration is about+2M(M-1)2+M(M-1),when there are K nonzero elements,the total number of the multiplication is

where[∗]denotes roundup.

According to[11],when there are K nonzero elements in source signals,the number of the multiplication in the OMP algorithm is

Thus the ratio of the OMP algorithm to the PP algorithm in computational complexity is about

According to[14],the number of multiplications in the CMP algorithm is

We can get the approximate ratio of complementary matching pursuit(CMP)to the PP algorithm in computational complexity:

The ASMP algorithm consists of inner iterations and outer iterations,and the backtracking refinement in ASMP would increase the number of iterations,so it usually costs more computational time than OMP,CMP and PP algorithms.Table 1 shows the value of γ1,γ2with different K,M,N.It can be seen that in most cases γ1and γ2are larger than 1 and on rare occasions γ2is a little less than 1,and in the UBSS problem,the M,N,K are usually not very large.Thus the computational complexity of the PP algorithm is lower than that of OMP and CMP algorithms in most cases.It can also be seen that γ1and γ2increase as K increases,which shows the PP algorithm will save more computational time when the source signal is insufficiently sparse.

Table 1 Value of γ1,γ2with different K,M,N

6.Simulation results

In this section,we compare the PP algorithm with OMP,BP and ASMP algorithms with the respective of the computational complexity and recovery accuracy.The simulations are divided into two categories.Simulation experiment 1 is concerned with the random signals,and we analyze the impact of the increase of SNR and the variation of sparsity on performance.Simulation experiment 2 evaluates the performance of algorithms for linear frequency modulation(LFM)signals.In the experiments,we use correlation coefficient[15]to measure recovery accuracy as follows:

where T is the total sampling time.s(t) =[s1(t),s2(t),...,sN(t)]Tdenotes the sourcesignalat time t and(t)=[1(t),2(t),...,N(t)]Tis the recovery signal at time t.The larger correlation coefficient is,the more accurate the algorithm will be.And the maximum value of γcoef(s,)is 1.In order to compare the computational complexity of different algorithms intuitively,we use the running time in the same simulation platform of different algorithms to measure computational complexity.Simulation results are the mean results of 2 000 experiments.

(i)Simulation experiment 1

The source signal is the random sparse signal with the sparsity p.The larger p is,the sparser the source signals are.There are three sensors and four sources and the mixing matrix is a random 3×4 matrix.Noise is modeled as a zero mean Gaussian noise.Fig.1 and Table 2 demonstrate the correlation coefficients and running time of four algorithms as a function of SNR when p=0.8 respectively.

Fig.1 Correlation coefficients with different SNRs when p=0.8

Table 2 Running time of the four algorithms with increase of SNR when p is 0.8 s

Fig.2 demonstrates the correlation coefficients of four algorithms versus sparsity p when SNR=30 dB.

Fig.2 Correlation coefficients with different p when SNR=0.8

In Fig.1,we observe that when SNR is less than 10 dB,the correlation coefficients of all the algorithms are nearly the same.When SNR is larger than 10 dB,the correlation coefficient of PP is larger than those of the other three algorithms.Fig.2 shows that when p＞0.4,the correlation coefficient of the PP algorithm is greater than that of other three algorithms.

The computational time of PP is reduced by 98%,86%over BP and ASMP algorithms respectively, and almost the same with the BP algorithm as shown in Table 2.Therefore,the PP algorithm can achieve a better balance between recovery accuracy and computational complexity than BP,OMP and ASMP algorithms.

(ii)Simulation experiment 2

The source signals are four LFM signals whose mathematical models aresi(t)=Aiexp(j(2π(fi-)t+)),where t is the sampling point,the sampling frequency is 400 MHz,Aiis amplitude,Bi=20 MHz,i= 1,2,3,4,and the number of the received signals is 3.f1=40 MHz,f2=40 MHz,f3=70 MHz,f4=70 MHz.The observation time is 2×10-6s,the mixing matrix is a 3×4 random matrix and noise is modeled as a zero mean Gaussian noise.It is easy to see that source signals are sparse in frequency domain and in most cases the number of nonzero sampling points is 2.Fig.3 and Table 3 demonstrate the correlation coefficients and the running time of four algorithms as a function of SNR.

Fig.3 Correlation coefficients with different SNRs

Table 3 Running time of the four algorithms with SNR s

It can be seen from Table 3 that the proposed PP algorithm reduces the computational time by 98%and 84%compared with BP and ASMP algorithms.It can be seen that,when SNR is smaller than 5 dB,the correlation coeffi cient of the PP algorithm is higher than that of the OMP algorithm,and is comparable with those of the BP algorithm and the ASMP algorithm.When SNR is larger than 5 dB,the correlation coefficient of the proposed PP algorithm is larger than those of the other three algorithms.

7.Conclusions

In this paper,a new sparse signal recovery algorithm,called the PP algorithm,is proposed.As compared with OMP,BP and ASMP algorithms,the PP algorithm achieves higher recovery accuracy,and has lower computational complexity.

[1]WANG X,HUANG Z,ZHOU Y.Underdetermined DOA estimation and blind separation of non-disjoint sources in time frequency domain based on sparse representation method.IEEE Systems Engineering and Electronics,2014:25(1):17–25.

[2]ZHANG L,YANG J,LU K,et al.Modi fi ed subspace method based on convex model for underdetermined blind speech separation.IEEE Trans.on Consumer Electronics,2014,60(2):225–232.

[3]HESSE C,JAMES C.On semi-blind source separation using spatial constraints with applications in EEG analysis.IEEE Trans.on Biomedical Engineering,2006,52(12):2525–2534.

[4]SHA Z,HUANG Z,ZHOU Y.Frequency-hopping signals sorting based on underdetermined blindsource separation.IET Communications,2013,7(14):1456–1464.

[5]ZHANG W,LIU J,SUN J.A new two-stage approach to underdetermined blindsource separation using sparse representation.IEEE International Conference on Acoustics,Speech and Signal Processing,2007:953–956.

[6]BAO G,YE Z,XU X.A compressed sensing approach to blind source separation of speech mixture based on a two-layer sparsity model.IEEE Trans.on Audio,Speech and Language Processing,2013,21(5):899–906.

[7]MALLEY M,NOVAK R.Near-optimal adaptive compressed sensing.IEEE Trans.on Information Theory,2014,60(7):4001–4012.

[8]CHEN S S,DONOHO D L,SAUNDERS M A.Atomic decomposition by basis pursuit.SIAM Review,2001,43(1):129–159.

[9]LIU E,TEMLYAKOV H.The orthogonal super greedy algorithm and application in compressed sensing.IEEE Trans.on Information Theory,2012,58(4):2040–2047.

[10]TROP J,GILBERT A.Signal recovery from random measurements via orthogonal matching pursuit.IEEE Trans.on Information Theory,2007,53(12):4655–4666.

[11]NEEDELL D,TROPP J.CoSaMP:iterative signal recovery from incomplete and inaccurate samples.Applied and Computational Harmonic Analysis,2009,26(3):301–321.

[12]DAI W,MILENKOVIC O.Subspace pursuit for compressive sensing signal reconstruction.IEEE Trans.on Information Theory,2009,55(5):2230–2249.

[13]WU H,WANG S.Adaptive sparsity matching pursuit algorithm for sparse reconstruction.IEEE Signal Processing Letters,2012,19(8):471–474.

[14]GAGANR,CHRISTINE G.Complementary matching pursuit algorithms for sparse approximation.Proc.of the 16th European Signal Processing Conference,2008:25–29.

[15]HE Z,XIE S,DING S.Convolutive blind source separation in the frequency domain based on sparse representation.IEEE Trans.on Audio,Speech and Language Processing,2007,15(5):1551–1563.