Abstract:
An affine projection or like algorithm is used to process a non-stationary signal. The affine projection algorithm creates an inverse matrix and includes a factor γ to avoid numerical instability in the inverse matrix. The factor γ is adaptively adjusted according to the characteristics of the signal.

Description:
FIELD OF THE INVENTION  
       [0001]     This invention relates to digital signal processing, and in particular to a method of processing a non-stationary signal, such as a speech signal, for example, for use in echo cancellation in speakerphone technology.  
       BACKGROUND OF THE INVENTION  
       [0002]     A typical echo cancellation application is shown in  FIG. 1 . The transmitted speech or signal is fed back to the received end becoming echo to the incoming speech. The echo reflection is caused by either impedance mismatch as in a network application or room reflection in an acoustical application. The echo will distort the incoming signal and speech.  
         [0003]     In telephony application, the echo is very annoying to the people on the phone. An adaptive algorithm is always applied to monitor the echo path and its variation so that it can be replicated by a linear filter and then subtracted from the incoming signal. For the telephony application, the speech signal is non-stationary and has a huge dynamic range.  
         [0004]     One adaptive algorithm is known as the Affine Projection Algorithm(APA), K. Ozeki and T. Umeda, “An Adaptive Filtering Algorithm Using an Orthogonal Projection to an Affine Subspace and Its Properties”, Electronics and Communications in Japan, Vol.67-A, No.5, 1984. This is a generalization of the well known normalized least mean square (NLMS) adaptive algorithm. It has a much faster convergence rate than the LMS algorithm and its fast implementation has similar computational complexity as the LMS algorithm.  
         [0005]     The LMS algorithm can be seen as a special case of APA where the projection order is 1. With the increase of projection order, the APA starts showing a dramatic increase of convergence speed. When the projection order is over 5, the increase of converge speed shows only a slight improvement. When the projection order increases to the order of filter size, the APA becomes the Recursive Least Square Algorithm (RLS).The APA is a very attractive algorithm which has the same complexity as LMS but is able to achieve RLS performance. In general, the convergence speed doubles with a doubling of projection order.  
         [0006]     The original APA needs more data memory than the NLMS algorithm and its complexity is equivalent to the N times of NLMS algorithm (here N is the order of projection. However, because of information and data reuse, many computation procedures can be saved and reused. This is the fundamental idea of the Fast Affine Projection (FAP) algorithm. However, the APA involves an Nth-order matrix inverse. There are different methods available for the matrix inverse which differentiates various algorithms with different performance, complexity and stability. The most important issue is that no matter what methods are, the problem of non-stability with large matrix inverse remains, especially with its fast operation and finite precision implementation. Singular value decomposition, J. Cioffi and T. Kailath, “Windowed Fast Transversal Filters Adaptive Algorithms with Normalization”, IEEE Trans. on ASSP, Vol:ASSP-33, No.3, June 1985, is the most stable implementation of matrix inverse and its only problem is its computational complexity.  
         [0007]     When dealing with a large matrix, there is always the singularity problem. With finite precision and fixed point implementation, the inverse of singular matrix will make the whole system very unstable because the smallest eigenvalue dominates the whole inverse matrix. In either APA or FAP algorithms, an identity matrix multiplied by a small constant is added to the inverse matrix so that the matrix will never become a singularity. However, this constant will cause performance degradation for the adaptive system. If this constant is too big, the performance degradation is severe and if this constant is too small, system stability cannot be maintained unless more data precision is employed in the echo cancellation.  
       SUMMARY OF THE INVENTION  
       [0008]     According to the present invention there is provided a method of processing a non-stationary signal comprising applying an adaptive echo cancellation algorithm that creates an inverse matrix to said non-stationary signal, and adaptively adjusting the condition of said matrix according to the characteristics of said signal.  
         [0009]     The adaptive echo cancellation algorithm may be an affine projection algorithm including a factor γ to avoid numerical instability in said inverse matrix. In this case the factor γ is adaptively adjusted.  
         [0010]     This method in accordance with the invention makes the algorithm much more stable and more able to track a speech signal with non-stationary properties. It also stabilizes the algorithm for higher order implementations.  
         [0011]     The novel method can be applied to other adaptive echo cancellation algorithm where the matrix inverse is required, such as the PNLMS (Proportionate Normalized Least-Mean-Squares) algorithms. The method can also be applied to an FAP algorithm with a forgetting factor α close to 1.  
         [0012]     In the method in accordance with the invention the APA can be used adaptively track the echo path change and then cancel the echo from the incoming signal.  
         [0013]     The novel method of the invention can, for example, be implemented in a digital signal processor.  
         [0014]     Test results show that with an adaptive matrix condition, the echo cancellation is similar to that using the singular value decomposition for low order APA and it is much better than constant approach. For a high order APA, the adaptive matrix condition approach outperforms singular value decomposition. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0015]     The invention will now be described in more detail, by way of example only, with reference to the accompanying drawings, in which:  
         [0016]      FIG. 1  shows a prior art echo cancellation circuit;  
         [0017]      FIG. 2  is a diagram showing a practical implementation of an echo canceller in accordance with a first embodiment of the invention;  
         [0018]      FIG. 3  is a simplified version of the implementation shown in  FIG. 2 ;  
         [0019]      FIG. 4  is a diagram showing a practical implementation of an echo canceller in accordance with a second embodiment of the invention;  
         [0020]      FIG. 5  shows a simulation echo speech signal and the echo path;  
         [0021]      FIG. 6  shows the convergence comparison between different algorithms;  
         [0022]      FIG. 7  shows an example with high order APA;  
         [0023]      FIG. 8  shows the simulated CSS (Composite Source Signal) echo with its echo path;  
         [0024]      FIG. 9  shows the echo cancellation performance comparison for the projection dimension N=4; and  
         [0025]      FIG. 9  shows the echo cancellation performance comparison for the projection dimension N=50. 
     
    
     DETAILED DESCRIPTION OF THE INVENTION  
       [0026]     The prior art Affine Projection Algorithm (APA) can be summarized as follows:  
         [0000]     Initialization: 
 
   x   (0)=   w   (0)=[0, 0, . . . , 0] T  
 
 for k&gt;0, 
 
   e   ( k )=   d   ( k )−   X     T ( k )   w   ( k ) 
 
   w   ( k +1)=   w   ( k )+μ   X   ( k )(   X     T ( k )   X   ( k )+γ   I   ) −1     e   ( k ) 
 
where 
 
   d   ( k )=[ d ( k ), d ( k −1), . . . ,  d ( k −N )] T  
 
 is the desired output. 
 
   w   ( k )=[ w   0 ( k ),  w   1 ( k ), . . . , w L-1 ( k )] T  
 
 is the weighting vector. 
 
   X   ( k )=[   x   ( k ),   x   ( k −1),. . . ,   x   ( k−N )]
 
   x   ( k )=[ x ( k ), x ( k −1), . . . , x ( k−L )] T  
 
 and 
 
 with x(k) being the reference signal at instant k. 0≦μ2 is the convergence step size, playing the same control role as in the LMS algorithm. γ is a small constant to avoid numerical instability in the matrix inverse. Normally, γ is chosen to be 30 dB below the peak value of the diagonal element of  X   T (k) X (k), or γ=max(x(k)*x(k))*L/1000. For the non-stationary signal with large magnitude variation such as speech signal, the choice of γ affects the adaptation convergence rate and final residual error. 
 
         [0027]     The detailed FAP algorithm is described in S. L. Gay and S. Tavathia, “The Fast Affine Projection Algorithm”, ICASSP&#39;95 Proceedings, pp.3023-3026; H. Ding, “A stable fast affine projection adaptation algorithm suitable for low-cost processors”, ICASSP 2000, Turkey, pp.360-363. In the derivation, one assumption has been used in the approximated fast residual echo vector calculation:  
           e   _     ⁡     (   k   )       ≈     [           e   ⁡     (   k   )                   (     1   -   μ     )     ⁢         e   _     _     ⁡     (     k   -   1     )               ]         
 
 (where  {overscore (e)} (k−1) is the upper N-1 elements of vector  e (k−1) This approximation is valid only if following condition satisfies γ&lt;&lt;λ i,k− 1, where λ i,k− 1 is the ith eigenvector of matrix    X     T ( k −1) X ( k −1). In general, if y is too big, the FAP performance will be affected and if γ is too small, the matrix inverse may become unstable, especial with finite precision implementation. 
 
         [0028]     Since the choice of γ affects the APA and FAP performance and stability, in accordance with the principles of the invention the value of γ is adaptively adjusted based on signal characteristics. Two embodiments are described.  
       FIRST EMBODIMENT  
       [0029]     The value of γ is set at 30 dB below the average power of last N signal set. In other words, γ is set at 30 dB below the average diagonal element of matrix  X   T (k) X (k). The equation can be written as  
       Υ   =       1     1000   ⁢   N       ⁢     ∑     diag   ⁡     (           X   _     T     ⁡     (   k   )       ⁢       X   _     ⁡     (   k   )         )               
 
         [0030]     This will ensure that the maximum eigen spread will be 1000N. To avoid the becoming zero or close to zero, in which case, the eigenspread is undefined (0/0), the above equation can be modified as:  
       Υ   =     max   ⁡     (         1     1000   ⁢   N       ⁢     ∑     diag   ⁡     (           X   _     T     ⁡     (   k   )       ⁢       X   _     ⁡     (   k   )         )           ,   c     )           
 
 where c is a constant representing the minimum precision in the signal representation. For the fixed point representation we can choose c=1, or a little larger than 1. With the definition of  X   T (k), the main components of the above equation can be written as:  
         ∑     diag   (           X   _     T     ⁡     (   k   )       ⁢       X   _     ⁡     (   k   )         )       =       ∑     i   =     k   -   N       k     ⁢     R   ⁡     (     0   ,   i     )             
     where     
         R   ⁡     (     0   ,   k     )       =       ∑     i   =   0     L     ⁢       x   ⁡     (     k   -   i     )       ⁢     x   ⁡     (     k   -   i     )               
 
         [0031]     The relation between R(0, k) and R(0, k−1) is 
 
 R (0 ,k )= R (0, k −1)+ x ( k ) x ( k )− x ( k−L− 1) x ( k−L− 1) 
 
 and the iterative relation is 
 
Σdiag(   X     T ( k )   X   ( k ))=Σdiag(   X     T ( k −1)   X (   k −1))+ R (0 ,k )− R (0 ,k −N −1) 
 
         [0032]     This will give an implementation diagram as shown in  FIG. 2 . Input signal x(k) is multiplied by itself in multiplier  10  and passed through delay unit  12  to adder  14  and from there to the input of delay unit  20  and adder  30 .  
         [0033]     The output of delay unit  20  is fed added to the output of adder  10  in adder  22  and fed to adder  14 .  
         [0034]     The output of delay unit  20  is passed to delay unit  24  and added in adder  26  to the output of delay unit  28 , which receives its input from adder  30 .  
         [0035]     The output of adder  30  is also fed to multiplier  32 , which receives as its other input 1/1000 N and then fed to maximization unit  18 . The other input of maximization unit  18  is a constant value c. The output of unit  18  is the adaptively adjusted factor γ.  
         [0036]     In most cases, L&gt;&gt;N and the diagonal elements in the matrix  X   T (k) X (k) would not have large variations from samples to samples. Therefore, as a simplified version, γ can be chosen as  
       Υ   =     max   ⁡     (         1     1000   ⁢   N       ⁢         ∑     i   =   0       L     ⁢       x   ⁡     (     k   -   i     )       ⁢     x   ⁡     (     k   -   i     )             ,   c     )           
 
 which gives a much simplified implementation as shown in  FIG. 3 . 
 
       SECOND EMBODIMENT  
       [0037]     This embodiment uses the estimated signal power R(0) as the basis and γ is chosen as 30 dB below that value. The estimated value R(0) is  
         R   ⁡     (   0   )       =       α   ⁢           ⁢     R   ⁡     (   0   )         +       (     1   -   α     )     ⁢       ∑     i   =   0     L     ⁢       x   ⁡     (     k   -   i     )       ⁢     x   ⁡     (     k   -   i     )                   
 
 where 0≦α≦1 is the forgetting factor. γ is calculated as:  
       Υ   =     max   ⁡     (         1   1000     ⁢     R   ⁡     (   0   )         ,       1   4000     ⁢         ∑     i   =   0       L     ⁢       x   ⁡     (     k   -   i     )       ⁢     x   ⁡     (     k   -   i     )             ,   c     )           
 
 where c is defined in the same way as before and the second term is used to control the initial value when x(k) just starts to become meaningful input. 
 
         [0038]     The practical implementation for this embodiment is shown in  FIG. 4 . The output of adder  14  is fed to subtractor  40 , having one input receiving the output of delay unit  42 .  
         [0039]     The output of adder  40  is fed to multiplier  44 , whose output is fed to the input of adder  46 .  
         [0040]     The output of adder  14  is also fed to the input of multiplier  50  whose output is fed to an input of maximization unit  18 . A third input of maximization unit  18  receives an output of multiplier  48 , which receives one input from the output of adder  46  and another input which represents a constant value.  
         [0041]      FIG. 5  shows the simulation echo speech signal and the echo path and  FIG. 6  shows the convergence comparison between different algorithms with N=4 for APA based algorithms.  
         [0042]     The second embodiment is applied for the modified APA algorithm (the simulations show that both the first and second embodiments perform very similarly). From  FIG. 6 , it can be seen that the APA algorithm barely outperforms the LMS algorithm and the modified APA performance has a much faster convergence rate and better ERLE. The modified APA based on the adaptive matrix condition estimation a similar echo cancellation performance as an SVD based APA. But the SVD implementation is much more complicated and very expensive for implementation.  
         [0043]      FIG. 7  shows an example with high order APA. Once again, this shows that the normal APA algorithm barely outperforms the LMS algorithm and the modified APA has a much faster convergence rate and better ERLE than the others.  
         [0044]      FIG. 8  shows the simulated CSS (Composite Source Signal) echo with its echo path. The CSS signal is a pseudo speech signal used to test echo cancellation performance in the ITU-T standard.  
         [0045]      FIGS. 9 and 10  show the echo cancellation performance comparison for the projection dimension N=4 and 50 respectively. Due to the fact that the CSS signal has almost constant magnitude, the normal APA performances is much better than the LMS algorithm. For N=4, the performance all three APA based algorithms is quite similar. This tells us that the APA with constant matrix condition value γ performs very well for a signal with stable magnitude, and the modified APA based on the adaptive matrix condition estimation outperforms other algorithms with high order projection.  
         [0046]     The schemes proposed in this paper can be applied to other adaptive echo cancellation algorithms where the matrix inverse is required, such as the PNLMS (Proportionate Normalized Least-Mean-Squares) algorithms described in S. Gay, “An Efficient, Fast Converging Adaptive Filter for Network Echo Cancellation”, Proc. Asilomar Conference on Signals, Systems and Computers, Pacific Grove, Calif., 1998; T. Gansler, S. L. Gay, M. M. Sondhi and J. Benesty, “Double-Talk Robust Fast Converging Algorithms for Network Echo Cancellation”, IEEE Trans. on Speech and Audio Processing, Vol.8, No.6, November 2000.pp.656-663; D. L. Duttweiler, “Proportionate Normalized Least-Mean-Squares Adaptation in Echo Cancelers”, IEEE Trans. on Speech and Audio Processing, Vol8, No.5, September 2000. pp.508-518.  
         [0047]     When the invention is applied to an FAP algorithm, similar performance improvement can be achieved. However, the forgetting factor α should be closer to 1 so that γ will only vary slowly. This is because the FAP uses an iterative method to create the matrix inverse and any dramatic change of the matrix with γ being varied will affect the stability of matrix inverse.  
         [0048]     All afore-mentioned references are herein incorporated by reference.