Abstract:
The present invention provides a method for blind separation of independent source signals. The preceding is accomplished by disclosing new techniques for dealing with the problem of blind separation of independent source signals. The method of the present invention consists of identifying a simple constrained criterion which stems directly from derived conditions for source separation. Using this criterion, observed and output signals are used to control a filtering matrix. After being initialized, the filtering matrix is updated and projected to the closest unitary matrix. Until the filtering matrix has converged, it continues to be updated and projected to the closest unitary matrix. Once convergence has occurred, the source signals have been recovered.

Description:
FIELD OF THE INVENTION 
     The present invention relates generally to the recovery of independent user signals simultaneously transmitted through a linear mixing channel, and more particularly, to a method and system for blind recovery of a number of independent user signals. 
     BACKGROUND OF THE INVENTION 
     In many signal processing applications, the sample signals provided by the sensors are mixtures of many unknown sources. The “separation of sources” problem is to extract the original unknown signals from these known mixtures. Generally, the signal sources as well as their mixture characteristics are unknown. Without knowledge of the signal sources other than the general statistical assumption of source independence, this signal processing problem is known in the art as the “blind source separation problem”. The separation is “blind” because nothing is known about the values of the independent source signals and nothing is known about the mixing process (which is assumed to be linear). 
     The blind separation problem is encountered in many familiar forms. For instance, the well-known “cocktail party” problem refers to a situation where the unknown (source) signals are sounds generated in a room and the known (sensor) signals are the outputs of several microphones. Each of the source signals is delayed and attenuated in some (time varying) manner during transmission from source to microphone, where it is then mixed with other independently delayed and attenuated source signals, including multipath versions of itself (reverberation), which are delayed versions arriving from different directions. A person, however, generally wishes to listen to a particular set of sound source while filtering out other interfering sources, including multi-path signals. 
     This signal processing problem arises in many contexts other than the simple situation where each of two mixtures of two speaking voices reaches one of two microphones. Other examples involving many sources and many receivers include the separation of radio or radar signals sensed by an array of antennas, sonar array signal processing, image deconvolution, radio astronomy, and signal decoding in cellular telecommunication systems. Those skilled in the signal processing arts have been eager to solve blind source separation problems because of their broad application to many communication fields. Solutions to date, however, are time intensive, require extensive computing power, and the source signal separation is not ideal. 
     SUMMARY OF THE INVENTION 
     The present invention broadly contemplates systems and methods for blind source separation based on the identification of a set of conditions which are necessary and sufficient for the separation of the source signals and which any method used for blind source separation must satisfy. Preferably, an optimization technique is used to enforce these conditions to separate and recover the source signals. This approach achieves better signal separation than known processes, and the optimization reduces the time and computing needed to separate the source signals. 
     In accordance with the present invention, a simple constrained criterion which stems directly from derived conditions for source separation is identified. Using this criterion, observed, i.e., first signals, and output signals are used to control a filtering matrix. After being initialized, the filtering matrix is updated and projected to the closest unitary matrix. Until the filtering matrix has converged, it continues to be updated and projected to the closest unitary matrix. Once convergence has occurred, the source signals have been recovered. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 shows a block diagram of the elements in a system utilized for blind separation of independent source signals in accordance with the present invention. 
     FIG. 2 shows a flow chart describing a method for blind separation of independent source signals in accordance with an illustrative embodiment of the present invention. 
     FIG. 3 shows a flow chart describing the preferred manner in which the projection to the closest unitary matrix is accomplished in accordance with an illustrative embodiment of the present invention. 
    
    
     DETAILED DESCRIPTION 
     Referring to the block diagram of the elements in a system utilized for blind separation of independent source signals in accordance with the present invention shown in FIG. 1, source or transmitted signals  10  are preferably p i.i.d. and mutually independent zero-mean discrete-time sequences a i (k) i=1, . . . ,p that share the same statistical properties. The source signals pass through transmission channel  20 , typically a p×q multiple-input-multiple-output (MIMO) linear memoryless channel which introduces interuser interference (IUI). Received signals  30  are inputted into pre-whitener  40 , resulting in observed or sensed signals  50 . Received signals  30  and sensed signals  50  may be written as Y(k) and Y′(k), respectively. Sensed signals  50  are then subsequently filtered by a q×p filtering matrix  60  whose outputs, output signals  70 , z j (k), j=1, . . . ,p should ideally match the source signals a i (k) when the proper filtering matrix  60  is used. To determine the proper filtering matrix  60 , output signals  70  are inputted, together with sensed signals  50 , into signal processor  80 , which controls filtering matrix  60 . 
     Source signal separation is achieved for both multiple user equalization and single user equalization when the filtering matrix  60  is selected such that it satisfies certain, in each case, conditions. Determining the proper equalizer setting in the case of single user equalization is a special case of multiple user equalization. The conditions to be satisfied in accordance with the present invention are as follows: 
     
       
         | K ( z   j ( k ))|=| K   a   |, j =1 , . . . , p   (C1) 
       
     
     
       
           E|z   j ( k )| 2 =σ a   2   , j =1 , . . . , p   (C2) 
       
     
     
       
           E ( z   i ( k ) z   j *( k ))=0 , i≠j   (C3) 
       
     
     E denotes the statistical expectation. For a detailed discussion of the derivation of these conditions, the reader is referred to the details which follow. 
     The only provision for an appendix is under rule 37 CFR 1.96 for reference to computer program listings which must be provided on CD. The reference to an Appendix was removed to avoid confusion and to avoid mis-treatment of pages 12-16 to ensure that they are properly printed as part of the specification.              {           max   G             F        (   G   )       =       ∑     j   =   1     p                          K        (     z   j     )                          subject                 to        :                 G   H        G     =   I                   (   1   )                                
     where I is the p×p identity matrix and H denotes the Hermitian (or else complex conjugate) transpose. The constraint in (1) comes from the fact that, according to (C2) and (C3), 
     
       
           E ( zz   H )=σ a   2   I   (2) 
       
     
     (It is assumed for simplicity that σ a   2 =1) In accordance with the present invention, (1) is the preferred maximization criterion. The global optima of this preferred criterion achieve blind recovery. 
     The p×q channel matrix is denoted by C and the q×1 channel output vector by Y(k). Although denoted Y(k), it is preferred this channel output be “pre-whitened” in both “space” and time by pre-whitener  40 . Thus, channel output Y(k) corresponds to sensed signal  50 . The sensed signal model is then: 
     
       
           Y ( k )= C   T   A ( k )+ n ( k )  (3) 
       
     
     where Y(k)=[y 1 (k) . . . y q (k)] T  and n(k) is the q×1 vector of additive noise samples. The receiver output, sensed signal  50 , can then be written as 
     
       
           z ( k )= W   T ( k ) Y ( k )= W   T ( k ) C   T   A ( k )+ n ′( k )= G   T ( k ) A ( k )+ n ′( k )  (4) 
       
     
     where W(k) and G(k)=CW(k) are the q×p receiver matrix and p×p global response matrix, respectively, and n′(k)=W T  (k)n(k) is the colored noise at the receiver output, all at time instant k. 
     FIG. 2 illustrates a flow chart of the steps executed by the signal processor  80  in determining the proper filtering matrix  60 , in accordance with the technique of the present invention, which may be referred to as a multi-user kurtosis (MUK) algorithm for blind source separation. Step S 20  involves initializing filtering matrix  60 , that is for k=0, setting W(0)=W 0 . This is accomplished by first computing the gradient of F(G) with respect to W. By writing F(G) as                F        (   G   )       =                  ∑     j   =   1     p                       sign        (     K        (     z   j     )       )            K        (     z   j     )                       =                  sign        (     K   a     )              ∑     j   =   1     p                     (       E               z   j          4       -     2        E   2                 z   j          2       -            E        (     z   j   2     )            2       )                                      
     and assuming symmetrical inputs (E(a i   2 (k)=0)), the gradient of F(G) with respect to W equals                ∇     (     F        (   G   )       )       =     4          ∑     j   =   1     p                     E        (                z   j          (   k   )            2            z   j          (   k   )              Y   *          (   k   )         )                   (   5   )                                
     At Step S 21 , for k&gt;0, it is determined if the matrix has converged. If the matrix has converged, the recovered signals outputted in Step S 22  are very close to the input signals. If convergence has not occurred, the filtering matrix is updated in Step S 23 . In accordance with the present invention, W(k) is updated in the direction of the in stantaneous gradient (dropping the expectation operator in (5)) as follows: 
     
       
           W ′( k +1)= W ( k )+μ sign( K   a ) Y *( k ) Z ( k )  (1) 
       
     
     where 
     
       
           Z ( k )=[| z   1 ( k )| 2   z   1 ( k ) . . . | z   p ( k )| 2   z   p ( k )]  (2) 
       
     
     and μ is the algorithm&#39;s step size. The orthogonality constraint needs to be satisfied at the next iteration of the algorithm: 
     
       
           G   H ( k +1) G ( k +1)= I   (3) 
       
     
     For this to be feasible, it is necessary for the channel output to be “pre-whitened” in both “space” and time. This corresponds to assuming the channel matrix C to be unitary, which is achieved by the pre-whitener  40  in FIG.  1 . Assuming that C is unitary, to satisfy the constraint (8) it suffices to satisfy 
     
       
           W   H ( k +1) W ( k +1)= I   (4) 
       
     
     Because there is no guarantee that W′(k+1) will satisfy the constraint (9), it must be hence transformed to a unitary W(k+1)=f(W′(k+1)). 
     Preferably, W(k+1) is chosen to be a q×p matrix which is as close as possible to W′(k+1) in the Euclidean sense. Dropping the time index for convenience, this can be achieved with an iterative procedure which satisfies the following criterion successively for j=1, . . . , p:              {           min     w   j               Δ        (     W   j     )       =                      W   j       -       W   j   ′                    2                     subject                 to        :                   W   l   H          W   j       =     δ   ij       ,     l   =   1     ,   …              ,   j                   (   10   )                                
     where ∥X∥ 2 =X H X is the squared Euclidean norm of vector X and where W′ is defined from 
     
       
           W′=[W′   1   . . . W′   p ]  (11) 
       
     
     The problem (10) can also be written as              {           min     w   j               Δ        (     W   j     )       =         (       W   j     -     W   j   ′       )     H          (       W   j     -     W   j   ′       )                   subject                 to        :                   W   l   H          W   j       =     δ   ij       ,     l   =   1     ,   …              ,   j                   (   12   )                                
     To solve the problem (12) the Lagrangian of Δ(W j ) is constructed, which equals                           Δ          (       W   j     ,     λ   j     ,   μ   ,   v     )       =                  Δ        (     W   j     )       -       λ   j          (         W   j   H          W   j       -   1     )       -                                  ∑     l   =   1       j   -   1                         μ   lj          Re        (       W   l   H          W   j       )           -       ∑     l   =   1       j   -   1                         v   lj          Im        (       W   l   H          W   j       )                           (   13   )                                
     where λ j , μ lj , v lj  are real scalar parameters. Setting the gradient of £ Δ (W j , λ j , μ, v) with respect to W j  to zero          (         ∂        Δ         ∂     W   j   *         =   0     )     ,                          
     the following are obtained for each j:                  W   j     -     W   j   ′     -       λ   j          W   j       -       ∑     l   =   1       j   -   1                         β   lj          W   l           =   0           (   14   )                                
     where β lj =½(μ lj +{square root over (−1)}v lj ). From (14) 
     
       
         
           
             
               
                 
                   { 
                   
                     
                       
                         
                           
                             λ 
                             j 
                           
                           = 
                           
                             1 
                             - 
                             
                               
                                 W 
                                 j 
                                 H 
                               
                                
                               
                                 W 
                                 j 
                                 ′ 
                               
                             
                           
                         
                       
                     
                     
                       
                         
                           
                             β 
                             ij 
                           
                           = 
                           
                             
                               - 
                               
                                 W 
                                 l 
                                 H 
                               
                             
                              
                             
                               W 
                               j 
                               ′ 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
                 
         
             
         
      
     
     is obtained, which gives                  (     1   -     λ   j       )          W   j       =       W   j   ′     -       ∑     l   =   1       j   -   1                         (       W   l   H          W   j   ′       )          W   l                   (   16   )                                
     According to (16),                W   j     ∝     (       W   j   ′     -       ∑     l   =   1       j   -   1                         (       W   l   H          W   j   ′       )          W   l           )             (   17   )                                
     where ∝ denotes “proportional”. In light of (17), and keeping in mind that W j   H W j =1, W j  (bringing back the time index) must be chosen as:                  W   j          (     k   +   1     )       =           W   j   ′          (     k   +   1     )       -       ∑     l   =   1       j   -   1                         (         W   l   H          (     k   +   1     )              W   j   ′          (     k   +   1     )         )            W   l          (     k   +   1     )                        W   j   ′          (     k   +   1     )       -       ∑     l   =   1       j   -   1                         (         W   l   H          (     k   +   1     )              W   j   ′          (     k   +   1     )         )            W   l          (     k   +   1     )                            (   18   )                                
     The updated filtering matrix is then projected to the closest unitary matrix in Step S 24 . FIG. 3 illustrates a flow chart of the steps taken in projecting the filtering matrix to the closet unitary matrix in Step S 24 . In Step S 30 , the norm of the first column of the updated filtering matrix  60  is corrected, that is            W   1          (     k   +   1     )       =         W   1   ′          (     k   +   1     )                     W   1   ′          (     k   +   1     )                                          
     In Step S 31  it is determined if there are more columns in the filtering matrix whose norm needs to be corrected. If not, a unitary filtering matrix is generated from the corrected columns in Step S 32 , according to equation (12). If so, the norm for the next column is corrected, that is: W j (k+1) is computed in accordance with Equation (18). 
     Steps S 30 - 33  correspond to a Gram-Schmidt orthogonalization of W′(k+1). The projection described by steps S 30 - 33  would reduce to the mere normalization of Step  530  in the case of a single user. Also, if the transmitted input is non-symmetrical, equation (6) should use instead of Z(k) in (7) the following vector: 
     
       
           Z ′( k )=[| z   1 ( k )| 2   z   1 ( k )−&lt; z   1   2 ( k )&gt; z   1 *( k ) . . . | z   p ( k )| 2   z   p ( k )−&lt; z   p   2 ( k )&gt; z   p *( k )]  (19) 
       
     
     where &lt; &gt; denotes empirical averaging. 
     The system of the present invention includes an input signal receiver, a filtering matrix, a signal processor, and a convergence detector, which may be implemented on a suitably programmed general purpose computer. These may also be implemented on an Integrated Circuit, part of an Integrated Circuit, or parts of different Integrated Circuits. Thus, it is to be understood that the invention may be implemented in hardware, software, or a combination of both. Accordingly, the present invention includes a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perform any of the method steps herein described for performing blind source separation of independent source signals. Such method steps could include, for example, receiving an input signal which is the result of an independent source signal having passed through a transmission channel; pre-whitening said input signal to obtain an observed signal; inputting said input signal into a filtering matrix to obtain a recovered signal; utilizing said observed signal and said recovered signal to update said filtering matrix; projecting said updated filtering matrix to the closest unitary matrix, outputting a signal representing said independent source signal when said recovered signal has converged, when said recovered signal has not converged, employing said closest unitary matrix as a new filtering matrix and continuing to input said input signal into a filtering matrix to obtain a recovered signal, utilizing said observed signal and said recovered signal to update said filtering matrix, projecting said updated filtering matrix to the closest unitary matrix, and outputting a signal representing said independent source signal when said recovered signal has converged. 
     Again, it is to be emphasized that any of the method steps, in any combination, can be encoded and tangibly embodied on a program storage device according to the present invention. 
     While only certain features of the invention have been illustrated and described herein, many modifications, substitutions, changes or equivalents will now occur to those skilled in the art. It is therefore intended to be understood that the appended claims are intended to cover all such modifications and changes that fall within the true spirit of the invention. 
     Multiple User Equalization 
     The equalizer output signals, z j (k), j=1, . . . , p can be written as a function of the source signals a i (k), i=1, . . . , p as follows:                    z   j          (   k   )       =         ∑     l   =   1     p                       g   jl            a   l          (   k   )           =       G   j   T          A        (   k   )             ,     j   =   1     ,   …              ,   p           (   20   )                                
     where A(k)=[a 1 (k) . . . a p (k)] T . G j  is the channel/equalizer cascade that contains the contribution of all the p channel inputs to the j-th equalizer output, and T denotes matrix transpose. Equation (20) can be also written in the familiar form 
     
       
           z ( k )=[ z   1 ( k ) . . .  z   p ( k )] T   =G   T   A ( k )  (21) 
       
     
     where              G   =       [       G   1                   …                   G   p       ]     =     [           g   11         ⋯         g   p1             ⋮       ⋮       ⋮             g     1      p           ⋯         g   pp           ]               (   22   )                                
     The above-described linear mixture setup can be used to model narrowband antenna-array systems, code division multiple access (CDMA) systems, or oversampled systems. 
     Each output&#39;s power and kurtosis may be expressed as:                  E        (            z                   (   k   )            2     )       =       σ   a   2            ∑     l   =   1     p                            g   jl          2           ,     j   =   1     ,   …              ,   p           (   23   )                   K        (     z   j     )       =       K   a            ∑     l   =   1     p                            g   jl          4           ,     j   =   1     ,   …              ,   p           (   24   )                                
     respectively K a  is the (unnormalized) kurtosis and σ a   2  the variance of each a j (k) (since they are all assumed to share the same distribution), defined as 
       K   a   =K ( a   j ( k )) 
     
       
         σ a   2   =E (| a   j ( k )| 2 )  (25) 
       
     
     where K(x)=E(|x| 4 )−2E 2 (|x| 2 )−|E(x 2 )| 2 . 
     Single User Equalization 
     Determining the proper equalizer setting in the case of single user equalization of intersymbol interference channels (ISI) channels is a special case of multiple user equalization. The single equalizer output may be written as z(k)=G T A(k), where G is a (possibly infinite-length) column vector containing the coefficients of the channel impulse response. Similar to equations (23) and (24), the output power and kurtosis for a single equalizer output are                E        (              z   j          (   k   )            2     )       =       σ   a   2            ∑     l   ∈   P                   g   l          2                 (   26   )                 K        (     z        (   k   )       )       =       K   a            ∑     l   ∈   P                   g   l          4                 (   27   )                                
     where P is a subset of the set of integers. Based on these expressions, the following set of conditions are necessary and sufficient for BE:              {                  K        (     z        (   k   )       )            =     K   a                   E               z        (   k   )            2       =     σ   a   2                     (   28   )                                
     which stems from the fact that                  ∑     l   ∈   P                   g   l          4       ≤       (       ∑     l   ∈   P                   g   l          2       )     2             (   29   )                                
     Equality in (29) holds only in the case that there exists a unique nonzero element g l  of unit magnitude for some lεP. Based on the condition (28), the following constrained optimization problem has been suggested for BE              {           max   G               F   SW          (   G   )       =          K        (     z   j     )                          subject                 to        :                              G   H        G     =   1                   (   30   )                                
     An important result about the behavior of criterion (11) is that the cost function F SW (G) is free of undesired local maxima on G H G=1 if G is infinite-length. This makes the corresponding algorithm for BE globally convergent to a “Dirac” solution, that is a solution with a single non-zero element (up to a scalar phase shift). Moreover, the popular Constant Modulus Algorithm (CMA) 2—2 can be seen as a special case of constrained optimization criteria when the input is sub-Gaussian, thus it reflects the global convergence property as well. 
     Linear Instantaneous Mixture Case 
     As is typical in blind equalization, the method of the present invention permits each transmitted signal to be recovered up to a unitary scalar rotation. This is due to the inherent inability of statistical blind techniques to distinguish between different rotated versions of the input signals. Therefore, blind recovery will be achieved if (after suitable reordering of the equalizer outputs) the following holds: 
     
       
           z   j ( k )= e   {square root over (−1)}φ     j     a   j ( k )  (31)  (31) 
       
     
     for some φ j ε[0,2π) and all jε{1, . . . , p}. 
     Based on (23) and (24), if each a i (k), i=1, . . . , p is an i.i.d. zero-mean sequence, {a i (k)}, {a j (k)} are statistically independent for i≠j and share the same statistical properties, then the following set of conditions are necessary and sufficient for the recovery of all the transmitted signals at the equalizer outputs: 
     
       
         | K ( z   j ( k ))|=| K   a   |, j =1 , . . . , p   (C1) 
       
     
     
       
           E|z   j ( k )| 2 =σ a   2   , j= 1 , . . . , p   (C2) 
       
     
     
       
           E ( z   i ( k ) z   j *( k ))=0 , i≠j   (C3) 
       
     
     To achieve perfect recovery (31) must hold, from which (C1) and (C2) follow immediately and (C3) follows since E(a i (k)a j *(k))=0 for i≠j (* denotes complex conjugate). 
     From (24) and (C1), the following is obtained:                  ∑     l   =   1     p                            g   jl          4       =   1           (   32   )                                
     From (23) and (C2), the following is obtained:                  ∑     l   =   1     p                            g   jl          2       =   1           (   33   )                                
     Therefore, equations (32) and (33) dictate, in the light of (29), that G j  must be of the form 
     
       
           G   j   =[. . . e   {square root over (−1)}φ     j   . . . ] T   (34) 
       
     
     where the single non-zero element can be at any position. Combining (C3) with (21) gives: 
     
       
           G   i   H   G   j =0 , i≠j   (35) 
       
     
     where H denotes the Hermitian (or else complex conjugate) transpose. According to (34) and (35), if G j &#39;s non-zero element is at position n, then G i  cannot have its non-zero element at position n. Hence, the p different “Dirac”-type vectors G j  of the form (34) will all contain their unique nonzero elements at different positions. This corresponds to the recovery of all the p different inputs a j (k). Therefore, after re-ordering, (31) is obtained and perfect signal separation (in the above sense) has been achieved.