Patent Publication Number: US-6343268-B1

Title: Estimator of independent sources from degenerate mixtures

Description:
BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The present invention relates to blind source separation problems and more particularly to estimating independent Auto Regressive (AR) processes from their sum. 
     2. Description of the Prior Art 
     A generic Blind Source Separation (BSS) problem is defined by the following: given m measurements x 1 , . . . ,x m  obtained from n independent signals (sources) s 1 , . . . ,s n , estimate the original signals through n estimators ŝ 1 , . . . , ŝ n , based on the time-series x 1 (·), . . . ,x m (·). 
     Current BSS literature addresses only the case when the number of sources is equal to the number of microphones, in example m=n. This is discussed by S. Amari in “Minimum Mutual Information Blind Separation”,  Neural Computation,  1996, by A. J. Bell and T. J. Sejnowski in “An Information-maximization Approach To Blind Separation And Blind Deconvolution”,  Neural Computation,  7:1129-1159, 1995, by J. F. Cardoso in “Infomax And Maximum Likelihood For Blind Source Separation”,  IEEE Signal Processing Letters,  4(4):112-114, April 1997, by P. Comon in “Independent Component Analysis, A New Concept?”,  Signal Processing,  36(3):287-314, 1994, by C. Jutten and J. Herault in “Blind Separation Of Sources, Part I: An Adaptive Algorithm Based On Neuromimetic Architecture”,  Signal Processing,  24(l):1-10, 1991, by B. A. Pearlmutter and L. C. Parra in “A Context-sensitive Generalization Of ICA”, In International Conference on  Neural Information Processing,  Hong Kong, 1996, by K. Torkkola in “Blind Separation Of Convolved Sources Based On Information Maximization”, In  IEEE Workshop on Neural Networks for Signal Processing,  Kyoto, Japan 1996. The case when the number of measurements m is strictly smaller than the number of sources n is called the degenerate case. 
     It is an object of the present invention to reconstruct independent signals from degenerate mixtures. More specifically, it is an object of the present invention to estimate independent Auto Regressive (AR) processes from their sum. 
     SUMMARY OF THE INVENTION 
     The present invention is a system that reconstructs independent signals from degenerate mixtures. More specifically, the present invention estimates independent Auto Regressive (AR) processes from their sum. The invention addresses the identification subsystem for such a degenerate case, particularly the case of two AR processes of known finite dimension (m=1 and n=2). 
     The present invention includes an identification system and an estimator. A mixture signal and noise are inputted into the system and through noise separation, a near pure signal is outputted. The identification system includes an ARMA identifier, a computation of autocovariance coefficients, an initializer and a gradient descent system. The estimator includes filtering. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 illustrates an example including the present invention. 
     FIG. 2 illustrates a block diagram of the present invention. 
     FIG. 3 illustrates a block diagram of the identification system. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The present invention is a system that reconstructs independent voices from degenerate mixtures by estimating independent Auto Regressive (AR) processes from their sum. The invention addresses the identification subsystem for such a degenerate case, particularly the case of two AR processes of known finite dimension (m=1 and n=2). 
     FIG. 1 illustrates an example utilizing the present invention. A mixture signal, such as two voices, is inputted into the system and through source identification and filtering by source separator  10 , an estimate of original sources are outputted from the present invention. The present invention estimates n unknown independent components S 1 , . . . ,s n  (sources) from observed signals x 1 , . . . ,x m  that are assumed to be linear combinations of the independent components. The present invention solves problems caused by degenerate source mixtures. 
     In the present invention the following assumptions are made: 
     1. There are two weakly stationary and independent sources S 1  and s 2  of zero mean, i.e. 
     
       
           E[s   j ( k ) s   j ( k−n )]= E[s   j ( l ) s   j ( l−n )], for every  k,l,j= 1,2 
       
     
     
       
           E[s   1 ( k ) s   2 ( l )]= E[s   1 ( k )]= E[s   2 ( l )]=0, for every  k,l   (1) 
       
     
     2. The source signals are obtained from two Auto Regressive (AR)(p) processes of order p as follows: 
     
       
           s   1 ( n )=− a   1   s   1 ( n− 1)− a   2   s   1 ( n− 2)− . . . − a   p   s   1 ( n−p )+ G   1   v   1 ( n ) 
       
     
     
       
           s   2 ( n )=− b   1   s   2 ( n− 1)− b   2   s   2 ( n− 1)− . . . − b   p   s   2 ( n−p )+ G   2   v   2 ( n )  (2) 
       
     
     where v 1  and v 2  are two independent unit variance white-noises, a 1 , . . . ,a p , G 1  and b 1 , . . . ,b p ,G 2  are the parameters of the first, respectively the second, AR(p) process. 
     3. The order p is assumed known. 
     4. The measurement is a scalar sum of the signals: 
     
       
           x ( n )= s   1 ( n )+ s   2 ( n )  (3) 
       
     
     and is known. 
     As stated above, the present invention estimates the original signals. The estimators ŝ 1 (·),ŝ 2 (·) that are constructed are optimal with respect to the conditional variance of the error. 
     The system solves the signal mixture separation problem. More specifically, a simultaneous record of two signals is given. The implementation of the system estimates each original signal. The system works best with unechoic voice signals, which can be well approximated by univariate AR processes of some order p. 
     As shown in FIG. 2, the major components of the present invention are the identifier  20  and the estimator  22 . There are two implementations of the present invention: an off-line implementation and an on-line implementation. Identifier  20  implements an identification system shown in FIG.  3  and estimator  22 , in one embodiment, is actually two Wiener filters. Identifier  20  produces an estimate of the unknown AR parameters â 1 , . . . ,â p ,Ĝ 1 ,{circumflex over (b)} 1 , . . . ,{circumflex over (b)} p ,Ĝ 2 . Estimator  22  gives the minimum variance estimates based on these parameters. 
     As shown in FIG. 3, sources s are linearly mixed in the environment. Identification system  32  receives input x from the mixing environment. Identification system  32  includes ARMA identifier  33 , initializer  34 , computation of autocovariance coefficients  35  and gradient descent system  36 . ARMA identifier  33 , initializer  34  and gradient descent system  36  receive parameters p, q, l and α from parameter block  37 . Outputted from identification system  32  are source estimates. 
     ARMA identifier  33  performs the operation described by step 2 below. Computation of autocovariance coefficients  35  performs the operation described by step 1 below. Initializer  34  performs the operation described by steps 3.1 through 3.4 and the gradient descent system  36  performs the operation of step 3.5 described below. 
     The off-line implementation assumes the entire measured signal is known from the start. Thus {x(n)} n  is known for n=1, . . . ,N and the method computes the estimates ŝ 1 (n)} n ,{ŝ 2 (n)} n  based on this knowledge. In the on-line implementation, at the moment t only the samples {x(n)} n  for n=1, . . . ,t are known. Based on this information, the system produces the estimates ŝ 1 (t+1),ŝ 2 (t+1) immediately after it receives x(t+1). Variations along this line are also considered. For instance the measured signal is split into time-windows of length M 0 ,            (         x   j          (   k   )       }     )         k   =   1     ,   …              ,     M   0         ,     j   =   1     ,   …              ,     N     M   0       ,     
              x   j          (   k   )       =     x        (         M   0        j     +   k     )         ,     k   =   1     ,   …              ,       M   0     .                     
     Then, for a given j 0 , immediately after the moment M 0 j 0 , the system produces the estimates 
     
       
           ŝ   1 (( j   0 −1) M   0   +k )} k=1, . . . ,M     0     ,{ŝ   2 (( j   0 −1) M   0   +k )} k=1, . . . ,M     0   . 
       
     
     Thus the length M 0  controls the time lag between the original signals and their estimates. 
     The present invention is based on the second order statistics of the measurement characterized through the sampled autocovariance coefficients computed in the off-line implementation through                  r   ^          (   l   )       =       1   N            ∑     k   =     l   +   1       N            x        (   k   )            x        (     k   -   l     )                     (   4   )                         
     In the on-line implementation, these coefficients are estimated through (24). 
     In the off-line implementation, identifier  20  performs the following steps with the following inputs: 
     
       
         { x ( n )} n=1, . . . ,N   ,P,  and constants:  L,q,α   0 , . . . ,α L ,ε 
       
     
     Step 1. Compute {circumflex over (r)}(l) for l=0, . . . ,L; 
     Step 2. Identify the measured signal as an ARMA(2p,p) process. 
     Thus two polynomials, one of degree p,                Q        (   z   )       =       ∑     l   =   0     p            c   l          z   l                 (   5   )                         
     and another of degree 2p,                P        (   z   )       =     1   +       ∑     l   =   1       2      p              d   l          z   l                   (   6   )                         
     are obtained. The zeros of P define the ARMA(2p,p) process poles. 
     Step 3. For each partition of the 2p poles identified at Step 2, into two admissible sets of p poles repeat: 
     3.1 Find P 1  and P 2  the p th  degree polynomials whose roots are, respectively, p poles from each partition. Set a 1   (0) , b 1   (0)  l=0, . . . ,p the coefficients of P 1 , P 2  as follows:                      P   1          (   z   )       =     1   +       ∑     l   =   1     p            a   l     (   0   )            z   l             ,       a   0     (   0   )       =   1            
                P   2          (   z   )       =     1   +       ∑     l   =   1     p            b   l     (   0   )            z   l             ,       b   0     (   0   )       =   1               (   7   )                         
     3.2 Solve the following linear system:                  [             ∑     l   =     -   p       p          h   l   2               ∑     l   =     -   p       p            h   l          g   l                     ∑     l   =     -   p       p            h   l          g   l                 ∑     l   =     -   p       p          g   l   2             ]     ·     [           x   1               x   2           ]       =     [             ∑     l   =     -   p       p            h   l          f   l                     ∑     l   =     -   p       p            g   l          f   l               ]             (   8   )                         
      where f 1 ,g 1 ,h 1  are obtained as coefficients of                    Q        (   z   )            Q        (     1   z     )         =       ∑     l   =     -   p       p            f   l          z   l                
                P   1          (   z   )              P   1          (     1   z     )         =       ∑     l   =     -   p       p            g   l          z   l                
                P   2          (   z   )              P   2          (     1   z     )         =       ∑     l   =     -   p       p            h   l          z   l                   (   9   )                         
     3.3 If x 1 &lt;0 or x 2 &lt;0 then this partition is not suitable and therefore move to the next partition. Go to Step 3. 
     3.4 Set 
     
       
           G   1   (0)   ={square root over (x 1 +L )},    G   2   (0)   ={square root over (x 2 +L )}   (10) 
       
     
     3.5 Using the gradient descent algorithm minimize the following criterion:                J        (       G   1     ,     G   2     ,   a   ,   b   ,   p   ,   L   ,   α   ,     r   ^     ,   q     )       =       ∑     l   =   0     L            α   1                   a   *   b   *       r   ^          (   l   )         -       ∑     k   =   l     p          (         G   1          b   k            Ψ   1          (     k   -   l     )         +       G   2          a   k            Ψ   2          (     k   -   l     )           )              q                 (   11   )                         
      where Ψ 1 (l),Ψ 2 (l), l=0, . . . ,p are obtained from the following triangular linear system,                      ∑     j   =   0     l            a   j            Ψ   1          (     l   -   j     )           =       G   1          δ     l   ,   0           ,     l   =   0     ,   …              ,   p          
                ∑     j   =   0     l            b   j            Ψ   2          (     l   -   j     )           =       G   2          δ     l   ,   0           ,     l   =   0     ,   …              ,   p             (   12   )                         
      with 
     
       
           a   0   =b   0 =1 
       
       
         
           
             
               
                 
                   
                     a 
                     * 
                     b 
                     * 
                     
                       
                         r 
                         ^ 
                       
                        
                       
                         ( 
                         l 
                         ) 
                       
                     
                   
                   = 
                   
                     
                       ∑ 
                       
                         k 
                         = 
                         0 
                       
                       
                         2 
                          
                         p 
                       
                     
                      
                     
                       
                         ∑ 
                         
                           j 
                           = 
                           0 
                         
                         k 
                       
                        
                       
                         
                           a 
                           j 
                         
                          
                         
                           b 
                           
                             k 
                             - 
                             j 
                           
                         
                          
                         
                           
                             r 
                             ^ 
                           
                            
                           
                             ( 
                             
                               l 
                               - 
                               k 
                             
                             ) 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
         
         
             
         
       
     
      α 1 , l=0, . . . ,L, are fixed positive weights and L,q are parameters. 
     Practically, the gradient descent step is performed as follows: 
     3.5.1 Set m=0 and compute 
     
       
           J   (old)   =J ( G   1   (m)   ,G   2   (m)   ,a   (m)   ,b   (m)   ,p,L,α,{circumflex over (r)},q )  (8) 
       
     
     3.5.2 Evaluate the gradient of J at G 1   (m) , G 2   (m) , a (m) , b (m)  in 
     
       
           Da   (m) =∇ a   J ( G   1   (m)   ,G   2   (m)   ,a   (m)   ,b   (m)   ,p,L,α,{circumflex over (r)},q ) 
       
     
     
       
           Db   (m) =∇ b   J ( G   1   (m)   ,G   2   (m)   ,a   (m)   ,b   (m)   ,p,L,α,{circumflex over (r)},q ) 
       
     
     
       
         
           
             
               
                 
                   
                     
                       DG 
                       1 
                       
                         ( 
                         m 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           ∂ 
                           J 
                         
                         
                           ∂ 
                           
                             G 
                             1 
                           
                         
                       
                        
                       
                         ( 
                         
                           
                             G 
                             1 
                             
                               ( 
                               m 
                               ) 
                             
                           
                           , 
                           
                             G 
                             2 
                             
                               ( 
                               m 
                               ) 
                             
                           
                           , 
                           
                             a 
                             
                               ( 
                               m 
                               ) 
                             
                           
                           , 
                           
                             b 
                             
                               ( 
                               m 
                               ) 
                             
                           
                           , 
                           p 
                           , 
                           L 
                           , 
                           α 
                           , 
                           
                             r 
                             ^ 
                           
                           , 
                           q 
                         
                         ) 
                       
                     
                   
                    
                   
                     
 
                   
                    
                   
                     
                       DG 
                       2 
                       
                         ( 
                         m 
                         ) 
                       
                     
                     = 
                     
                       
                         
                           ∂ 
                           J 
                         
                         
                           ∂ 
                           
                             G 
                             2 
                           
                         
                       
                        
                       
                         ( 
                         
                           
                             G 
                             1 
                             
                               ( 
                               m 
                               ) 
                             
                           
                           , 
                           
                             G 
                             2 
                             
                               ( 
                               m 
                               ) 
                             
                           
                           , 
                           
                             a 
                             
                               ( 
                               m 
                               ) 
                             
                           
                           , 
                           
                             b 
                             
                               ( 
                               m 
                               ) 
                             
                           
                           , 
                           p 
                           , 
                           L 
                           , 
                           α 
                           , 
                           
                             r 
                             ^ 
                           
                           , 
                           q 
                         
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
         
         
             
         
       
     
     3.5.3 Set t=1 and repeat 
       ax=a   (m)   −t·Da   (m)   
     
       
         
           bx=b 
           (m) 
           −t·Db 
           (m) 
         
       
     
     
       
         
           Gx 
           1 
           =G 
           1 
           (m) 
           −t·DG 
           1 
           (m) 
         
       
     
     
       
         
           Gx 
           2 
           =G 
           2 
           (m) 
           −t·DG 
           2 
           (m) 
         
       
     
     
       
           J   x   =J ( Gx   1   ,Gx   2   ,ax,bx,p,L,α,{circumflex over (r)},q ) 
       
     
     
       
           t=t/ 2  (16) 
       
     
      until J x &lt;J (old)    
     3.5.4 Set 
     
       
           a   (m+1)   =a   (m) −2 tDa   (m)   
       
     
     
       
           b   (m+1)   =b   (m) −2 tDb   (m)   
       
     
     
       
           G   1   (m+1)   =G   1   (m) −2 tDG   1   (m)   
       
     
     
       
           G   2   (m+1)   =G   2   (m) −2 tDG   2   (m)   
       
     
     
       
           J   (new)   =J   x   (17) 
       
     
     3.5.5 If (J (old) −J (new) )/J (old) ε set 
     
       
           m=m+ 1 
       
     
     
       
           J   (old)   =J   (new)   (18) 
       
     
      and go to 3.5.2. 
      Else set 
     
       
         
           a 
           (partition) 
           =a 
           (m+1) 
         
       
     
     
       
         
           b 
           (partition) 
           =b 
           (m+1) 
         
       
     
     
       
         
           G 
           1 
           (partition) 
           =G 
           1 
           (m+1) 
         
       
     
     
       
         
           G 
           2 
           (partition) 
           =G 
           2 
           (m+1) 
         
       
     
     
       
           J   opt (partition)= J   (new)   (19) 
       
     
     Step 4. Choose the partition that gives the smallest criterion, i.e.                  optim   —        partition     =     arg                     min   partition            J   opt          (   partition   )                   (   20   )                         
     Set: 
       â=a   option     —     partition   
     
       
         
           {circumflex over (b)}=b 
           option 
           
             — 
           
           partition 
         
       
     
     
       
         
           Ĝ 
           1 
           =G 
           1 
           option 
           
             — 
           
           partition 
         
       
     
     
       
           Ĝ   2   =G   2   option     —     partition   (21) 
       
     
     Outputs: â,{circumflex over (b)},Ĝ 1 ,Ĝ 2 . 
     This identification procedure uses a number of parameters. Experimentally, the following values work well: 
     
       
           L= 3 p   
       
     
     
       
           q= 2 
       
     
     
       
         α l =1,  l= 0 , . . . ,L   
       
     
     
       
         ε=10 −3   (22) 
       
     
     For the purposes of the Step 3, a set of complex numbers is said admissible if it is invariant under the complex conjugation operator, i.e. if z is in the set then {overscore (z)} should belong to the set too. 
     Once the identification procedure has been done, the next step is the linear filtering of the signal x performed in estimator  20  of FIG.  2 . The following will describe the off-line filtering block of estimator  20 . The present invention uses Wiener filters to estimate ŝ 1  and ŝ 2 . The Wiener filters are the minimum variance linear estimators and have been widely used in the literature as noise reduction filters. For the AR sources estimation, the filters are defined through the following relations:                    H   1          (   z   )       =         G   1   2              P   ^     2          (   z   )                P   ^     2          (     1   z     )                 G   ^     1   2              P   ^     2          (   z   )                P   ^     2          (     1   z     )         +         G   ^     2   2              P   ^     2          (   z   )                P   ^     2          (     1   z     )                    
              H   2          (   z   )       =         G   2   2              P   ^     1          (   z   )                P   ^     1          (     1   z     )                 G   ^     1   2              P   ^     2          (   z   )                P   ^     2          (     1   z     )         +         G   ^     2   2              P   ^     1          (   z   )                P   ^     1          (     1   z     )                       (   23   )                         
     where                P   ^     1          (   z   )       =     1   +       ∑     l   =   1     p              a   ^     l          z   l             ,           P   ^     2          (   z   )       =     1   +       ∑     l   =   1     p              b   ^     l            z   l     .                             
     Since (23) defines non-causal filters, the algorithm uses a forward-backward implementation for each filter, based on the spectral factorization theorem. 
     The following will describe the on-line implementation and the on-line identification algorithm. For the on-line problem, the present invention uses an adapting mechanism of the identification procedure. Suppose the estimates â (m) ,{circumflex over (b)} (m) ,Ĝ 1   (m) ,Ĝ 2   (m)  of the parameters a,b,G 1  and G 2  based on the first m samples are given and the m+1 st  sample x(m+1) is known. The on-line algorithm describes how these estimates are adapted. 
     Input: â (m) ,{circumflex over (b)} (m) ,Ĝ 1   (m) ,Ĝ 2   (m) ,{circumflex over (r)} (m) (l) and x(m+1) 
     Step 1. Update the autocovariance coefficients using an exponential-window, 
     
       
           {circumflex over (r)}   (m+1) ( l )=(1− w ) {circumflex over (r)}   (m) ( l )+ wx ( m+ 1) x ( m+ 1− l ),  l= 0 , . . . ,L   (24) 
       
     
     Step 2. Perform one gradient descent step 3.5 from the off-line algorithm using the updated values of the autocovariance coefficients:                t   *     =     arg                     min     t   =     2     -   n                J        (         G   1     (   m   )       -     t   ·     DG   1         ,       G   2     (   m   )       -     t   ·     DG   2         ,       a     (   m   )       -     t   ·   Da       ,       B     (   m   )       -     t   ·   Db       ,   p   ,   L   ,   α   ,       r   ^       (     m   +   1     )       ,   q     )                   (   25   )                         
     Update the new estimates: 
       a   (m+1)   =a   (m)   −t*Da   
     
       
         
           b 
           (m+1) 
           =b 
           (m) 
           −t*Db 
         
       
     
     
       
         
           G 
           1 
           (m+1) 
           =G 
           1 
           (m) 
           −t*DG 
           1 
         
       
     
     
       
           G   2   (m+1)   =G   2   (m)   −t*DG   2   (26) 
       
     
     Step 3. Set m=m+1 and go to Step 1. 
     Outputs: â (m+1) ,{circumflex over (b)} (m+1) ,Ĝ 1   (m+1) ,Ĝ 2   (m+1) ,{circumflex over (r)} (m+1) (l) 
     The parameter w at Step 1 represents the adapting rate of the algorithm. Usually it is taken very close to 0. 
     The following will describe on-line filtering. The filtering block adapts its parameters (through (23) at every M 0  samples. Depending on the application, either a forward-backward implementation is used (if the processes are slow), or the forward filtering part only is used (if the processes are fast). In the latter case, the forward filter is the causal (i.e. stable) part of (23). 
     The present invention solves problems caused by degenerate source mixtures. The present invention can be utilized in any application requiring speech enhancement and noise reduction such as hearing aids, mobile phones, stationary phones and speech recognition and detection systems. The present invention also has applications with source separation in EKG, EEG, magneto-encephalograms and fetal monitoring. 
     It is not intended that this invention be limited to the hardware or software arrangement or operational procedures shown disclosed. This invention includes all of the alterations and variations thereto as encompassed within the scope of the claims as follows.