WO2000038319A1 - Stable adaptive filter and method - Google Patents

Stable adaptive filter and method Download PDF

Info

Publication number
WO2000038319A1
WO2000038319A1 PCT/CA1999/001068 CA9901068W WO0038319A1 WO 2000038319 A1 WO2000038319 A1 WO 2000038319A1 CA 9901068 W CA9901068 W CA 9901068W WO 0038319 A1 WO0038319 A1 WO 0038319A1
Authority
WO
WIPO (PCT)
Prior art keywords
adaptive filter
step size
linear equations
solving
calculator
Prior art date
Application number
PCT/CA1999/001068
Other languages
English (en)
French (fr)
Inventor
Heping Ding
Original Assignee
Nortel Networks Corporation
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Priority claimed from US09/218,428 external-priority patent/US6754340B1/en
Application filed by Nortel Networks Corporation filed Critical Nortel Networks Corporation
Priority to JP2000590294A priority Critical patent/JP2002533970A/ja
Priority to CA002318929A priority patent/CA2318929A1/en
Priority to EP99973514A priority patent/EP1057259A1/en
Publication of WO2000038319A1 publication Critical patent/WO2000038319A1/en

Links

Classifications

    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03HIMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00Adaptive networks
    • H03H21/0012Digital adaptive filters
    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03HIMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00Adaptive networks
    • H03H21/0012Digital adaptive filters
    • H03H21/0043Adaptive algorithms
    • HELECTRICITY
    • H03ELECTRONIC CIRCUITRY
    • H03HIMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H21/00Adaptive networks
    • H03H21/0012Digital adaptive filters
    • H03H21/0043Adaptive algorithms
    • H03H2021/0049Recursive least squares algorithm
    • H03H2021/0052Recursive least squares algorithm combined with stochastic gradient algorithm
    • H03H2021/0054Affine projection

Definitions

  • the present invention relates to adaptive filters, and in particular, to fast affine projection
  • FAP FAP adaptive filters providing a stability of operation, and methods of stable FAP adaptive filtering.
  • Adaptive filtering is a digital signal processing technique that has been widely used in technical areas such as, e.g., echo cancellation, noise cancellation, channel equalization, system identification and in products like, e.g., network echo cancellers, acoustic echo cancellers for full-duplex handsfree telephones and audio conference systems, active noise control, data communications systems.
  • an adaptive filter The characteristics of an adaptive filter are determined by its adaptation algorithm.
  • the choice of the adaptation algorithm in a specific adaptive filtering system directly affects the performance of the system.
  • N MS normalized least mean square
  • LMS least mean square
  • the NLMS algorithm converges slowly with colored training signals like the speech, an important class of signals most frequently encountered in many applications such as telecommunications.
  • the performance of systems incorporating NLMS adaptive filters very often suffers from the slow convergence nature of the algorithm.
  • Other known algorithms proposed so far are either too complicated to implement on a commercially available low- cost digital signal processor (DSP) or suffer from numerical problems.
  • DSP digital signal processor
  • FAP fast affine projection
  • a method of adaptive filtering comprising the steps of:
  • (c) updating the filter coefficients comprising: determining auto-correlation matrix coefficients from a reference input signal, and solving at least one system of linear equations whose coefficients are the auto-correlation matrix coefficients, the system being solved by using a descending iterative method having an inherent stability of its operation, the results of the solution being used for updating the filter coefficients and the number of systems of linear equations to be solved being dependent on the normalized step size; (d) repeating the steps (b) and (c) required number of times.
  • the normalized step size may be chosen to be equal to any value from 0 to 1 depending on the application. In the majority of applications, it is often set to be close to unity or equal to unity. Conveniently, the normalized step size is within a range from about 0.9 to 1.0. Another convenient possibility is to set the normalized step size within a range from about 0.7 to 1.0.
  • the step of solving at least one system of linear equations comprises solving one system of linear equations only. Alternatively, in some applications, e.g., when one needs to keep misadjustment low after convergence, it is required to set the normalized step size substantially less than unity, e.g.
  • the step of solving at least one system of linear equations comprises solving N systems of linear equations, with N being a projection order.
  • N being a projection order.
  • a problem of finding the inverse of an auto-correlation matrix which is inherent for other known methods is reduced to a problem of solving a system of linear equations based on the autocorrelation matrix.
  • the system is solved by one of descending iterative methods which provide inherent stability of operation due to an intrinsic feedback adjustment. As a result inevitable numerical errors are not accumulated.
  • a steepest descent and conjugate gradient methods are used respectively to determine the first column of the inverse auto-correlation matrix, taking into account that the normalized step size is close to unity.
  • a steepest descent or conjugate gradient method is used to determine coefficients of the inverse auto-correlation matrix by recursively solving N systems of linear equations having decrementing orders. It corresponds to the case of the normalized step size being not close to unity.
  • the forth embodiment of the invention avoids determining the inverse of the auto-correlation matrix. Instead, a system of linear equations is solved by using a conjugate gradient method resulting in a solution that can be used directly to determine an updating part of the filter coefficients .
  • other known descending methods e.g. steepest descent, Newton's method, PARTAN, quasi-Newton' s method or other known iterative descending methods may also be used. Conveniently, the steps of the method may be performed by operating with real value or complex value numbers .
  • the method described above is suitable for a variety of applications, e.g. echo cancellation, noise cancellation, channel equalization, system identification which are widely used in products such as network echo cancellers, acoustic echo cancellers for full-duplex handsfree telephones and audio conference systems, active noise control systems, data communication systems.
  • an adaptive filter comprising: a filter characterized by adaptive filter coefficients; means for updating the filter coefficients, including means for setting a normalized step size, the updating means comprising: a correlator for determining auto-correlation matrix coefficients from a reference input signal, and a calculator for solving at least one system of linear equations whose coefficients are the auto-correlation matrix coefficients, the system being solved by using a descending iterative method having an inherent stability of its operation, the results of the solution being used for updating the filter coefficients and the number of systems of linear equations to be solved being dependent on the normalized step size.
  • the calculator is an iterative calculator.
  • the calculator is a steepest descent or a conjugate gradient calculator.
  • it may be a calculator performing a Newton's or quasi- Newton's method, a PARTAN calculator, or another known iterative descending calculator providing an inherent stability of operation.
  • the filter and the updating means are capable of operating with real numbers. Alternatively, they may be capable of operating with complex numbers.
  • the normalized step size may be chosen to be equal to any value from 0 to 1 depending on the application. In the majority of applications, the adaptive filter is often set with the normalized step size close to unity or equal to unity. Conveniently, the normalized step size is within a range from about 0.9 to 1.0. Another convenient possibility is to set the normalized step size within a range from about 0.7 to 1.0. For the normalized step size close to unity, the calculator provides iterative solution of one system of linear equations only at each time interval. Alternatively, in some applications, e.g., when one needs to keep misadjustment after convergence low, it is required to set the normalized step size substantially less than unity, e.g. less than about 0.7.
  • the calculator provides solutions of N systems of linear equations, with N being a projection order.
  • determining of the inverse auto-correlation matrix may be performed by solving N systems of linear equations having decrementing orders .
  • the adaptive filter as described above may be used for echo cancellation, noise cancellation, channel equalization, system identification or other applications where adaptive filtering is required.
  • the adaptive filter and method described above have an advantage over known FAP adaptive filters by providing a stability of operation.
  • the problem caused by error accumulation in matrix inversion process existing in known FAP filters is solved in the present invention by using iterative descending methods.
  • the matrix inversion operation is reduced to a solution of a corresponding system of linear equations based on the auto-correlation matrix.
  • the iterative descending methods, used for the solution of the above system provide an inherent stability of operation due to an intrinsic feedback adjustment. As a result, inevitable numerical errors are not accumulated, thus providing stability of adaptive filtering.
  • Figure 1 is a block diagram of an adaptive echo cancellation system
  • FIG. 2 is a block diagram of an adaptive filter according to the first embodiment of the invention.
  • Figure 3 is a block diagram of a steepest descent calculator embedded in the filter of Fig. 2 ;
  • Figure 4 is a block diagram of a conjugate gradient calculator embedded in an adaptive filter according to a second embodiment of the invention;
  • Figure 5 is a block diagram of an adaptive filter according to a third embodiment of the invention
  • Figure 6 is a flow-chart illustrating an operation of a steepest descent calculator embedded in the adaptive filter of Fig. 5;
  • Figure 7 is a flow-chart illustrating an operation of a conjugate gradient calculator embedded in the adaptive filter of Fig. 5;
  • Figure 8 is a block diagram of an adaptive filter according to a fourth embodiment of the invention.
  • Figure 9 is a block diagram of a conjugate gradient calculator embedded in the adaptive filter of Fig. 8.
  • d(n) and X(n) stand for column vectors, and bold-faced ones, like X(n), are matrices.
  • d(n) stands for an N-1 vector consisting of the N-1 upper most elements of the N vector d(n)
  • d(n) stands for an N-1 vector consisting of the N-1 lower most elements of the N vector d(n) .
  • a superscript "T” stands for the transposition of a matrix or vector.
  • FIG. 1 presents a block diagram of an adaptive echo cancellation system 10 with an embedded adaptive filter 100, the echo cancellation being chosen as an exemplary representation of a wide class of adaptive filtering applications .
  • a digitally sampled far-end reference input signal x(n) is supplied to the adaptive filter 100 and to an echo path 14 producing an unwanted signal u(n), the signal being an echo of x(n) through the echo path 14.
  • the echo path 14 can be either a long electrical path, e.g. in a telecommunication network, or an acoustical path, e.g. in a room.
  • An echo canceller may be used together with a telecomminication network switch or a speaker phone.
  • the unwanted signal u(n) is mixed up with the wanted near-end signal s(n) in a summer 16 to produce a response signal d(n) .
  • the response signal d(n) is sent to another summer 18 together with an echo estimate signal y(n) generated by the adaptive filter 100.
  • the summer 18 subtracts y(n) from d(n) 4
  • the adaptive filter Since the echo path is constantly changing, the adaptive filter must be able to continuously adapt to the new echo path. Therefore the goal is to produce the echo estimate signal y(n) as close to u(n) as possible, so that the latter is largely cancelled by the former, and e(n) best resembles s(n) .
  • the output signal e(n), called the error- signal is then transmitted to the far-end and also sent to the adaptive filter 100 which uses it to adjust its coefficients .
  • far-end and “near-end” may need to be interchanged.
  • x(n) in Figure 1 is actually the near-end signal to be transmitted to the far-end
  • d(n) in Figure 1 is the signal received from the telephone loop connected to the far-end.
  • the terminology used above is based on the assumption that x(n) is the far-end signal and d(n) is the signal perceived at the near-end, it is done solely for convenience and does not prevent the invention from being applied to other adaptive filter applications with alternate terminology.
  • d(n) and X(n) stand for column vectors, and bold-faced ones, like X(n), are matrices .
  • d(n) stands for an N-1 vector consisting of the N- 1 upper most elements of the N vector d(n)
  • d(n) stands for an N-1 vector consisting of the N-1 lower most elements of the N vector d(n) .
  • a superscript "T” stands for the transposition of a matrix or vector.
  • L-dimensional column vectors are defined as the reference input vector and the adaptive filter coefficient vector respectively, where L is the length of the adaptive filter:
  • Equation 1 The part for convolution and subtraction, which derives the output of the adaptive echo cancellation system, can then be expressed as
  • Equation 2 Equation 2 where the superscript "T” stands for transpose of a vector or matrix.
  • W(n + 1) W(n) + 2 ⁇ (n)e(n)X(n)
  • Equation 3 ⁇ (n) is called the adaptation step size, which controls the rate of change to the coefficients, is a normalized step size, and 6, being a small positive number, prevents ⁇ (n) from going too big when there is no or little reference signal x(n).
  • the computations required in the NLMS filter include 2L+2 multiply and accumulate (MAC) operations and 1 division per sampling interval.
  • LMS least mean square
  • the affine projection method is a generalization of the NLMS method. With N being a so-called projection order, we define
  • Equation 4 Equation 4 where d(n) and e(n) are N vectors and X(n) is an LxN matrix. Usually N is much less than L, so that X(n) having more a "portrait” rather than a “landscape” shape. Note that e(n) in Equation (4) is the a priori error vector; all its ele- ments, including e(n-l), ..., e(n-N+l), depend on W(n), as indicated in Equation (5) below.
  • Equation 5 W(n) is defined in Equation (1) .
  • W(n+1) W(n) + ⁇ X(n) ⁇ (n)
  • Equation 6 Equation 6 where I is the NxN identity matrix, and and ⁇ play similar roles as described with regards to Equation 3. is the normalized step size which may have a value from 0 to 1, and very often is assigned a unity value, ⁇ is a regular- ization factor that prevents R(n), the auto-correlation matrix, from becoming ill-conditioned or rank-deficient, in which case P(n) would have too big eigenvalues causing instability of the method. It can be seen that an NxN matrix inversion operation at each sampling interval is needed in the AP method.
  • the AP method offers a good convergence property, but computationally is very extensive. It needs 2LN+0(N 2 ) MACs at each sampling interval. For example, for N equal to 5, which is a reasonable choice for many practical applica- tions, the AP is more than 5 times as complex as the NLMS.
  • the FAP method consists of two parts:
  • Equation (7) (a) An approximation which is shown in Equation (7) below and certain simplifications to reduce the computational load.
  • the approximation in Equation (7) uses the scaled posteriori errors to replace the a priori ones in Equation
  • the matrix inversion may be performed by using different approaches.
  • One of them is a so-called “sliding windowed fast recursive least square (FRLS) " approach, outlined in US patent 5,428,562 to Gay, to recursively calculate the P(n) in Eq. 6. This results in a total requirement of computations to be 2L+14N MACs and 5 divi- sions.
  • the matrix inversion lemma is used twice to derive P(n) at sampling interval n, see, e.g. Q. G. Liu, B. Champagne, and K. C.
  • Equation (6) the first expression in Equation (6) shows that the coefficient vector W(n) will no longer be updated properly. That is, W(n) can be updated in wrong directions, causing the adaptive filtering system to fail.
  • a proposed remedy is to periodically re-start a new inversion process, either sliding windowed FRLS or conventional RLS based, in parallel with the old one, and to replace the old one so as to get rid of the accumulated numerical errors in the latter. While this can be a feasible solution for high-precision DSPs such as a floating-point processor, it is still not suitable for fixed-point DSP implementations because then the finite precision numerical errors would accumulate so fast that the re-starting period would have to be made unpractically short.
  • Equation 8 where P(n) is the very first, i.e., left most, column of the matrix P(n) .
  • P(n) is the very first, i.e., left most, column of the matrix P(n) .
  • P(n) is the very first, i.e., left most, column of the matrix P(n) .
  • P(n) is the very first, i.e., left most, column of the matrix P(n) .
  • P(n) is the very first, i.e., left most, column of the matrix P(n) .
  • Q.G. Liu cited above that, even with an ⁇ slightly less than that range, say about 0.7, the approximation is still acceptable.
  • N rather than all the N 2 , elements of P(n) .
  • Equation 10 Equation 10 where R(n) is symmetric and positive definite according to its definition Equation (9), and b is an Nvector with all its elements zero except the very first, which is unity.
  • determining of an updating part of the filter coefficients may be performed either by a direct solving for £(n) (second line of Eq.
  • a method of adaptive filtering implemented in an adaptive filter 100 according to the first embodiment of the invention includes an iterative "steepest descent” technique to iteratively solve the Equation (10) .
  • steepest descent is a technique that seeks the minimum point of a certain quadratic function iteratively. At each iteration (the same as sampling interval in our application) , it takes three steps consec- utively :
  • the steepest descent reaches the unique minimum of the quadratic function, where the gradient is zero, and continuously tracks the minimum if it moves. Details about the steepest descent method can be found, for example, in a book by David G.
  • the implied quadratic function is as follows ip T (n)R(n)P(n)-P T b
  • Equation 12 Equation 12
  • R(n) must be symmetric and positive definite in order for the steepest descent technique to be applicable, this happens to be our case. Seeking the minimum, where the gradient vanishes, is equivalent to solving Equation (10) .
  • the steepest descent is also able to track the minimum point if it moves, such as the case with a non-stationary input signal X(n) .
  • the stable FAP (SFAP) method which uses the steepest descent technique includes the following steps: Initialization:
  • R(n) R(n-l) + ⁇ (n) ⁇ T (n)- ⁇ (n-L) ⁇ T (n-L)
  • R(n) is the first column of R(n)
  • R(n) is an N-1 vector that consists of the N-1 lower most elements of the N vector R(n)
  • ⁇ (n) is an N-1 vector that consists of the N-1 upper most elements of the N vector ⁇ (n) .
  • Equation (15) determines P(n) based on P(n-l) and the new incoming data X(n) only, without examining how well a P actually approximates R _1 (n) . Therefore inevitable numerical errors will accumu- late and eventually make the system collapse.
  • the feedback provided by a stable descending method, used in our invention uses Equation (15) to examine how well P(n-l), or the needed part of it, approximates R " (n), or its corresponding part. Then the adjustments are performed in Equations (16) and (17) accordingly to derive P(n), or the needed part of it. As just mentioned, this examination is done by evaluating g(n) in Equation (15) as the feedback error.
  • Equation (15) is the gradient of the implied quadratic function (Equation (15) )
  • ⁇ (n) is the optimum step size for parameter vector adjustment, which is made in Equation (17) .
  • the total computational requirement of the Stable FAP method according to the first embodiment of the invention is 2L+2N 2 +7N-1 MACs and 1 division. Note, that for the steepest descent technique to work adequately for the purpose of adaptive filtering, the projection order N has to be chosen to assure that the steepest descent converges faster than the adaptive filter coefficients do. The required pre-determined value of N will depend on a particular adaptive filtering application.
  • An adaptive filter 100 according to the first embodi- ment of the invention and operating in accordance with the method described above is shown in Figure 2. It includes a filter 102 characterized by adaptive filter coefficients W(n), and means 104 for updating the coefficients, the means being set with a normalized step size ⁇ close to its maximal value, i.e. unity.
  • the filter 102 is a finite impulse response (FIR) filter which receives a Table 1
  • the updating means 104 includes a correlator 106 for recursively determining an auto-correlation signal presented in the form of auto-correlation matrix coefficients R(n) based on the reference input signal x(n), and a calculator 108 for generating projection coefficients P(n), the projection coefficients being part of the coefficients of the inverse of the auto-correlation matrix.
  • the calculator 108 defines projection coefficients by using an iterative steepest descent method having an inherent stability of operation as illustrated in detail above.
  • the projection coefficients are used within updating means 104 for generation the auxiliary filter adapta- 24 tion signal f(n) and an echo estimate correction signal EC(n) (see Equation (34) below) .
  • the latter is used together with the provisional echo estimate PR(n) to produce the echo estimate signal y(n) .
  • a convention in Fig. 2 is the use of a thick line to represent the propagation of a matrix or vector signal, i.e., with more than one component, and the use of a thin line to stand for a scalar signal propagation.
  • a correlator 106 determines the autocorrelation matrix R(n) in accordance with the Eq. 14 using the current and past x(n) samples.
  • An “lj.(n) calculator” 110 calculates H.(n) based on Eq. 22, and as shown in Fig. 2, T ⁇ (n) is not used by the updating means 104 until the next sampling interval.
  • the filter 102 produces the convolutional sum W (n)X(n) .
  • ⁇ N _ ⁇ (n-l) is obtained from by putting the latter through a unit delay element 111, providing a delay of one sampling interval, and further multiplied by the step size in a Multiplier 113. The result is used for updating the adaptive filter coefficients in (Eq. 18) . n.
  • ⁇ (n-l) is dot-multi- plied with part of R(n) by a Dot multiplier 112, and the result is further multiplied by a multiplier 114 with the step size ⁇ to form the correction term to be added to W T (n)X(n) by the summer 116 to form the filter output y(n) (Equation (19) ) .
  • the summer 18 calculates the error, or the output, e(n), as in Equation (20) .
  • the scalar-vector multiplier 118 derives £(n) in accordance with Equation (21) .
  • a steepest descent calculator 108 is shown in detail in Figure 3. Thick lines represent the propagation of a matrix or vector signal, i.e., with more than one component, and the use of a thin line stands for a scalar sig- nal propagation.
  • the autocorrelation matrix R(n) and the vector P(n-l) which is a part of the estimated inverse of R(n-l) are multiplied in a Matrix-vector multiplier 130.
  • the vector product is fur- ther subtracted by a constant vector [1 0 ... 0] ⁇ in a Summer 132 to produce the gradient vector g(n), which contains the feedback error information about using P(n-l) as the estimated inverse of R(n).
  • Equation (15) The squared norm of g(n) is then found by dot- multiplying g(n) with itself in a Dot multiplier 134. It is used as the numerator in calculating ⁇ (n) in Equation 16.
  • a Matrix-vector multiplier 136 finds the vector product between the autocorrelation matrix R(n) and the gradient vector g(n). This vector product is then dot-multiplied with g(n) in another Dot multiplier 138 to produce the denominator in calculating ⁇ (n) in Equation (16) .
  • This denominator is reciprocated in a Reciprocator 140, and then further scalar-multiplied with the aforementioned numerator in scalar multiplier 142 to produce ⁇ (n) . This is the only place where any division operation is performed. Finally, ⁇ (n) is multiplied with the gradient g(n) in a scalar-vector multiplier 144 to form the correction term to P(n-l). This correction term is then subtracted from P(n-l) in a Vector Summer 146 to derive P(n) in accordance with Equation (17). P(n-l) is obtained from P(n) by using a unit delay element 148, providing a delay of one sampling interval .
  • the first one is a floating point module
  • the second one is a 16-bit fixed-point DSP implementation.
  • a floating-point module simulating the NLMS acoustic echo canceller design in Venture, a successful full-duplex handsfree telephone terminal product by Nortel Networks Corporation, and a bench mark, floating- point module that repeats a prior art FAP scheme by Q .
  • G. Liu, B. Champagne, and K. C. Ho Bell-Northern Research and INRS-Telecommunications, Universite du Quebec), "On the Use of a Modified Fast Affine Projection Algorithm in Subbands for Acoustic Echo Cancellation," pp.
  • the source ones are speech files with Harvard sentences (Intermediate Refer- ence System filtered or not) sampled at 8 KHz and a white noise file. Out of the source files certain echo files have been produced by filtering the source ones with certain measured, 1200-tap, room impulse responses. These two sets of files act as x(n) and d(n) respectively.
  • the major simulation results are as follows.
  • the output e(n) in the steepest descent embodiment converges approximately at the same speed as the bench mark prior art FAP and reaches the same steady state echo cancellation depth as the prior art FAP and NLMS.
  • the SFAP according to the first embodiment of the invention outperforms NLMS filter; with speech training, it converges in about 1 second while it takes the NLMS filter about 7 to 8 seconds to do so.
  • a method of adaptive filtering according to a second embodiment of the present invention uses an iterative "conjugate gradient” technique to iteratively solve the Equation (10), the corresponding calculator being shown in Figure 4.
  • Conjugate gradient is a technique that also seeks the minimum point of a certain quadratic function iteratively. Conjugate gradient is closely related to the steepest descent scheme discussed above. It differs from the steepest decent in that it is guaranteed to reach the minimum in no more than N steps, with N being the order of the system. That is, conjugate gradient usually converges faster than the steepest descent. At each iteration (the same as sampling interval in out application) , the conjugate gradient takes five steps consecutively:
  • conjugate gradient modifies the negative gradient to determine an optimized direction.
  • the scheme reaches the unique minimum of the quadratic func- tion, where the gradient is zero, in no more than N steps.
  • the conjugate gradient technique also continuously tracks the minimum if it moves, such as the case with non-stationary input signal x(n) . Details about the conjugate gradient algorithm can be found, for example, in a book by David G. Luenberger (Stanford University), Linear and Non- linear Programming, Addison-Wesley Publishing Company, 1984.
  • Equation (11) whose gradient with respect to P(n) is also Equation (12) .
  • R(n) must be symmetric and positive definite in order for the conjugate gradient technique to apply, this happens to be our case. Seeking the minimum, where the gradient vanishes, is equivalent to solving Equation (10) .
  • the conjugate gradient is also able to track the minimum point if it moves, such as the case with non-stationary input signal X(n) .
  • the SFAP method according to the second embodiment which uses the conjugate gradient technique, includes the following steps: Initialization:
  • Equation 25 where (n) is defined in Equation (23) above, and determining projection coefficients by solving the system of linear Equations (10) using the conjugate technique, the projection coefficients being first column coefficients of the inverse of the auto-correlation matrix:
  • Equation 37 where R(n) is the first column of R(n), R(n) is an N-1 vector that consists of the N-1 lower most elements of the N vector R(n), and ⁇ (n) is an N-1 vector that consists of the N-1 upper most elements of the N vector ⁇ (n) .
  • Equations (26), (27), (28) , (31) and (32) respectively correspond to the five steps of the conjugate gradient technique discussed earlier in this section.
  • g( ) is the gradient of the implied quadratic function
  • ⁇ (n) is the optimum factor for updating the direction vector s(n).
  • ⁇ (n) is the optimum step size for parameter vector adjustment, which is made in Equation
  • the total computational requirement of the Stable FAP method according to the second embodiment of the invention is 2L+2N 2 +9N+1 MACs and 1 division. It should be also ensured that the conjugate gradient converges fast enough so that the adaptive filter coeffients converge.
  • An adaptive filter according to the second embodiment of the invention is similar to that of the first embodi- ment shown in Figure 2 except for the calculator 108 now operating in accordance with the conjugate gradient technique and being designated by numeral 208 in Figure 4.
  • the conjugate gradient calculator 208 embedded in the adaptive filter of the second embodiment is shown in detail in Figure 4. Thick lines represent the propagation of a matrix or vector signal, i.e., with more than one component, and the use of a thin line stands for a scalar signal propagation.
  • the autocorrelation matrix R(n) and the vector P(n-l), part of the esti- mated inverse of R(n-l) are multiplied in a Matrix-vector Multiplier 210.
  • the resulted vector product is subtracted by a constant vector [1 0 ...
  • a Summer 212 to produce the gradient vector g(n), which contains the feedback error information about using P(n-l) as the estimated inverse of R(n).
  • the Matrix-vector Multiplier 210 and the Summer 212 implement the Equation (26) above.
  • the gradient g(n) is further dot-multiplied with b(n-l), an auxiliary vector found in the last sampling interval, in a Dot Multiplier 214.
  • the resulted scalar product is multiplied by r srs (n-l) in a Multiplier 216, to produce ⁇ (n), a factor to be used in adjusting s(n-l), the direction vector for adjusting J
  • Equation (27) The part of the dia- gram described in this paragraph implements Equation (27) shown above. With ⁇ (n), g(n), and s(n-l) available, s(n-l) is then updated into s(n) by using yet another unit delay element 222, with a delay of one sampling interval, scalar- vector Multiplier 224 and Vector Summer 226 which imple- ment operations shown in Equation (28) above.
  • the auxiliary vector b(n), to be used in the next sampling interval is calculated as the product between R(n) and s(n) in another Matrix-vector Multiplier 230. This implements Equation (29) above.
  • the vector b(n) is then dot-multi- plied with s(n) in yet another Dot multiplier 232, and the scalar product is reciprocated in a Reciprocator 234, to produce r srs (n) (Equation (30)) . This is where the only division operation is.
  • g(n) and s(n) are dot-multiplied, and the result, being a scalar product, is multiplied with -r srs (n) to derive ⁇ (n), thus implementing Equation (31) above.
  • ⁇ (n) is available, it is multiplied with s(n) in another scalar-vector Multiplier 240 to form the correction term to P(n-l), which is then added to P(n-l) in a Vector Summer 242 in order to derive P(n) (Equation (32) above) .
  • the output e(n) in the conjugate gradient embodiment converges approximately at the same speed as the bench mark prior art FAP and reaches the same steady state echo cancellation depth as the bench mark prior art FAP and NLMS .
  • the SFAP according to the second embodiment of the invention also ourperformes NLMS filter in terms of convergence speed.
  • a method of adaptive filtering according to a third embodiment of the present invention provides adaptive filtering when the normalized step size has any value from 0 to 1. It updates the adaptive filter coefficients by iteratively solving a number of systems linear equations having decrementing orders to determine the inverse auto- correlation matrix in a manner described below.
  • Equation 39 Equation 39
  • Equation 40 This means that P is also the inverse of R. Since the inverse of a matrix is unique, the only possibility is
  • Equation 42 Equation (42) can be written in a scalar form
  • Equation 43 where r ik (n) is the element of R(n) on row i and column k, and p k :(n) the element of P(n) on row k and column j, and dy is defined as
  • Equation (45) coincides with Equation (10) derived earlier and applied to the first and second embodiments of the invention.
  • Equation 46 The right hand side of Equation (45) or Equation (46) tells that P(n) is the left-most column of P(n) and, based on Equation (41), P (n) is also the upper-most row of P(n) . According to the first and second embodiments of the invention discussed above, this part will cost "2N 2 +3N" MACs and 1 division with steepest descent or "2N +5N+2" MACs and 1 division with conjugate gradient.
  • Equation (47) can be re-arranged to become
  • these N-1 unknowns can be uniquely determined by only N-1 equations.
  • Equation 49 Equation (49) has the same format as Equation (45) except that the order is reduced by one. Equation (49) can also be solved by using either of the two approaches presented above, costing "2 (N-1) 2 +4 (N-1) MACs and 1 division with steepest descent” or "2 (N-1) 2 +6 (N-1) +2 MACs and 1 division with conjugate gradient," where the added "(N-1)” in each of the two expressions accounts for the extra computations needed to calculate the right hand side of Equation (49) .
  • Equation 49 Equation 49 Equation (49) has the same format as Equation (45) except that the order is reduced by one. Equation (49) can also be solved by using either of the two approaches presented above, costing "2 (N-1) 2 +4 (N-1) MACs and 1 division with steepest descent” or "2 (N-1) 2 +6 (N-1) +2 MACs and 1 division with conjugate gradient," where the added "(N-1)” in each of the two expressions accounts for
  • Equation (45) and Equation (49) are just special cases of Equation (50), and (p kj ( n ) , k-j , ' +l, ... , N- 1 ⁇ found in recursion step j form a column vector P;(n) , which consists of the lower N- elements of the j'th (0 ⁇ ⁇ N- 1) column of P(n) .
  • the process of Equation (50) will take ⁇ divisions and
  • Equation 51 for steepest descent method, and N divisions and [2N 2 + 5N + 2] + [2(N - l) 2 + 6(N - 1) + 2] + [2(N - 2) 2 + 7(N - 2) + 2] + ...
  • the SFAP method according to the third embodiment of the invention includes the following steps: Initialization:
  • Equation 54 Updating the adaptive filter coefficients in sampling interval n including the steps shown in Equation 55 below.
  • designations used in Equation (55), are as follows: (n) is defined in Equation (23) above, R(n) is the first column of R(n), R(n) is an N-1 vector that consists of the N-1 lower most elements of the N vector R(n), and ⁇ (n) is an N-1 vector that consists of the N-1 upper most elements of the Nvector ⁇ (n) .
  • any division operation in the 2nd expression of Equation (55) is not performed if the denominator is not greater than zero, in which case a zero is assigned to the quotient .
  • the filter 300 also differs from the filter 100 by the following features: the normilized step size may have any value from 0 to 1.0, the calculator 308 now has more extended structure for consec- utively determining columns of the inverse auto-correlation matrix in accordance with the steepest descent technique, and an e(n) calculator 320 is added.
  • the routine 402 sets an initial value to index j (block 404) which is submitted together with the auto-correlation matrix R(n) (block 406) to a projection coefficient column calculator (block 408) .
  • the calculator provides a steepest descent iteration in accordance with Equation (50) for the current value of index j, thus updating the corresponding column of projection coefficients from the previous sampling interval (block 408) .
  • the updated column of the projected coefficients is sent to a storage means (routine 410, block 412) to be stored until the other columns of P(n) are calculated.
  • thick lines represent the propagation of a matrix or vector signal, i.e., with more than one component, and the use of a thin line stands for a control propagation.
  • the steepest descent calculator 308 may be replaced with the conjugate calculator.
  • the corresponding structure is illustrated by a flow-chart 500 shown in Figure 7 where the blocks similar to that ones of Figure 6 are designated by same refer- ence numerals incremented by 100. It operates in a manner described above with regard to Figure 6.
  • a method of adaptive filtering according to a fourth embodiment of the present invention also provides adaptive filtering when the normalized step size has any value from 0 to 1. It updates the adaptive filter coefficients by iteratively solving a number of systems linear equations which avoid an explicit matrix inversion performed in the third embodiment of the invention. The details are described below.
  • the second equation from the set of Equations (6) which is reproduced for convenience in Equation (56) below, is equivalent to
  • Equation 56 it is possible to obtain £(n), required for updating the adaptive filter coefficients, directly from the set of linear Equations (56) , which are solved again by one of the descending iterative methods .
  • the SFAP method of the fourth embodiment of the invention includes the following steps:
  • T ⁇ -r srs (k+l)g ⁇ (k+l) Nx (N+l)
  • ⁇ t (k+l) ⁇ t (k) + ⁇ s(k+l) NxN
  • An adaptive filter 600 according to a fourth embodiment of the invention is shown in detail in Figure 8. It includes a filter 602 characterized by adaptive filter- coefficients W(n), and means 604 for updating the coefficients, the means being set with a normalized step size ⁇ having any value in a range from 0 to 1.0.
  • the filter 602 is a finite impulse response (FIR) filter which receives a reference input signal x(n) and an auxiliary signal f(n) used for updating the coefficients, and generates a provisional echo estimate signal PR(n).
  • FIR finite impulse response
  • the updating means 604 includes a correlator 606 for recursively determining an auto-correlation signal presented in the form of auto-correlation matrix coefficients R(n) based on the reference input signal x(n), an ⁇ (n) calculator 608 and an e(n) calculator 620 for corresponding calculation of vectors £_(n) and e(n) .
  • the calculator 608 defines ⁇ (n) by using an iterative conjugate gradient method having an inherent stability of operation as illustrated in detail above.
  • the projection coefficients are used within updating means 604 for generation the auxiliary filter adaptation signal f(n) and an echo estimate correction signal EC(n) . The latter is used together with the provisional echo estimate PR(n) to produce the echo estimate signal y(n).
  • Fig. 8 thick lines represent propagation of a matrix or vector signal, i.e., the signal with more than one component, and the use of a thin line stands for a scalar signal propagation.
  • a correlator 606 determines the autocorrelation matrix 4
  • R(n) in accordance with the first formula of Eq. (59) using the current and past x(n) samples.
  • An "n(n) calculator" 610 calculates n . (n) based the last formula of Eq. (59) , and as shown in Fig. 8, rj.(n) is not used by the updating means 104 until the next sampling interval.
  • the filter 602 produces the convolutional sum W (n)X(n) .
  • ⁇ N _ 1 (n-l) is obtained from ⁇ N -iC 11 ) ky putting the latter through a unit delay element 611, providing a delay of one sampling interval, and further multiplied by the step size ⁇ in a Multiplier 613.
  • the £(n) calculator 608 solves the sixth equation of Eq. (59) for £(n) by a conjugate gradient method, thus providing sufficient data for updating the adaptive filter coefficients (Eq. 6, first formula) .
  • Thick lines represent the propagation of a matrix or vector signal, i.e., with more 4
  • the calculator 608 additionally includes an output switch 754 which automatically opens- at the beginning of the sampling interval and closes at the end of N conjugate gradient iterations. Modifications described with regard to the first two embodiments are equally applicable to the third and fourth embodiments of the invention.
  • an adaptive filter and a method providing a stability of adaptive filtering based on feedback adjustment are provided.

Landscapes

  • Filters That Use Time-Delay Elements (AREA)
  • Cable Transmission Systems, Equalization Of Radio And Reduction Of Echo (AREA)
PCT/CA1999/001068 1998-12-22 1999-11-10 Stable adaptive filter and method WO2000038319A1 (en)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2000590294A JP2002533970A (ja) 1998-12-22 1999-11-10 安定適応フィルタおよびその方法
CA002318929A CA2318929A1 (en) 1998-12-22 1999-11-10 Stable adaptive filter and method
EP99973514A EP1057259A1 (en) 1998-12-22 1999-11-10 Stable adaptive filter and method

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
US09/218,428 US6754340B1 (en) 1998-12-22 1998-12-22 Stable adaptive filter and method
US09/218,428 1999-07-16
US09/356,041 1999-07-16
US09/356,041 US6788785B1 (en) 1998-12-22 1999-07-16 Stable adaptive filter and method

Publications (1)

Publication Number Publication Date
WO2000038319A1 true WO2000038319A1 (en) 2000-06-29

Family

ID=26912898

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CA1999/001068 WO2000038319A1 (en) 1998-12-22 1999-11-10 Stable adaptive filter and method

Country Status (4)

Country Link
EP (1) EP1057259A1 (ja)
JP (1) JP2002533970A (ja)
CA (1) CA2318929A1 (ja)
WO (1) WO2000038319A1 (ja)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2371191A (en) * 2001-01-11 2002-07-17 Mitel Corp Detecting double-talk and path changes in echo cancellers
DE10329055A1 (de) * 2003-06-27 2005-01-20 Infineon Technologies Ag Verfahren und Vorrichtung zur Echokompensation
WO2006017934A1 (en) * 2004-08-17 2006-02-23 National Research Council Of Canada Adaptive filtering using fast affine projection adaptation
CN106105032A (zh) * 2014-03-20 2016-11-09 华为技术有限公司 用于自适应滤波器的系统和方法
CN106209023A (zh) * 2016-07-28 2016-12-07 苏州大学 基于数据重用的非负自适应滤波方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0684706A1 (en) * 1993-11-05 1995-11-29 Ntt Mobile Communications Network Inc. Replica producing adaptive demodulating method and demodulator using the same
US5568411A (en) * 1994-07-25 1996-10-22 National Semiconductor Corporation Method and apparatus for using polarity-coincidence correlators in LMS adaptive filters
EP0751619A1 (en) * 1995-06-30 1997-01-02 Nec Corporation Noise cancelling method and noise canceller

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0684706A1 (en) * 1993-11-05 1995-11-29 Ntt Mobile Communications Network Inc. Replica producing adaptive demodulating method and demodulator using the same
US5568411A (en) * 1994-07-25 1996-10-22 National Semiconductor Corporation Method and apparatus for using polarity-coincidence correlators in LMS adaptive filters
EP0751619A1 (en) * 1995-06-30 1997-01-02 Nec Corporation Noise cancelling method and noise canceller

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
KAZUHIKO OZEKI ET AL: "AN ADAPTIVE FILTERING ALGORITHM USING AN ORTHOGONAL PROJECTION TO AN AFFINE SUBSPACE AND ITS PROPERTIES", ELECTRONICS AND COMMUNICATIONS IN JAPAN,US,SCRIPTA TECHNICA. NEW YORK, vol. 67-A, no. 5, 1 May 1984 (1984-05-01), pages 19 - 27, XP000563345 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2371191A (en) * 2001-01-11 2002-07-17 Mitel Corp Detecting double-talk and path changes in echo cancellers
GB2371191B (en) * 2001-01-11 2005-06-15 Mitel Corp Double-talk and path change detection using a matrix of correlation coefficients
DE10329055A1 (de) * 2003-06-27 2005-01-20 Infineon Technologies Ag Verfahren und Vorrichtung zur Echokompensation
DE10329055B4 (de) * 2003-06-27 2005-10-13 Infineon Technologies Ag Verfahren und Vorrichtung zur Echokompensation
WO2006017934A1 (en) * 2004-08-17 2006-02-23 National Research Council Of Canada Adaptive filtering using fast affine projection adaptation
US7436880B2 (en) 2004-08-17 2008-10-14 National Research Council Of Canada Adaptive filtering using fast affine projection adaptation
CN106105032A (zh) * 2014-03-20 2016-11-09 华为技术有限公司 用于自适应滤波器的系统和方法
EP3108579A4 (en) * 2014-03-20 2017-09-20 Huawei Technologies Co., Ltd. System and method for adaptive filter
CN106209023A (zh) * 2016-07-28 2016-12-07 苏州大学 基于数据重用的非负自适应滤波方法
CN106209023B (zh) * 2016-07-28 2018-06-26 苏州大学 基于数据重用的非负自适应滤波方法

Also Published As

Publication number Publication date
EP1057259A1 (en) 2000-12-06
CA2318929A1 (en) 2000-06-29
JP2002533970A (ja) 2002-10-08

Similar Documents

Publication Publication Date Title
US6788785B1 (en) Stable adaptive filter and method
Enzner et al. Frequency-domain adaptive Kalman filter for acoustic echo control in hands-free telephones
Comminiello et al. Nonlinear acoustic echo cancellation based on sparse functional link representations
US7171436B2 (en) Partitioned block frequency domain adaptive filter
Gil-Cacho et al. Nonlinear acoustic echo cancellation based on a sliding-window leaky kernel affine projection algorithm
Nascimento et al. Adaptive filters
Khan et al. Fractional LMS and NLMS algorithms for line echo cancellation
Huang et al. Practically efficient nonlinear acoustic echo cancellers using cascaded block RLS and FLMS adaptive filters
CN109102794A (zh) 基于凸组合的m估计成比例类仿射投影的回声消除方法
EP1782533A1 (en) Adaptive filtering using fast affine projection adaptation
Schrammen et al. Efficient nonlinear acoustic echo cancellation by dual-stage multi-channel Kalman filtering
Van Vaerenbergh et al. A split kernel adaptive filtering architecture for nonlinear acoustic echo cancellation
EP1314247B1 (en) Partitioned block frequency domain adaptive filter
WO2000038319A1 (en) Stable adaptive filter and method
Ramdane et al. Partial update simplified fast transversal filter algorithms for acoustic echo cancellation
Panda et al. A VSS sparseness controlled algorithm for feedback suppression in hearing aids
Okhassov et al. Cost-Effective Proportionate Affine Projection Algorithm with Variable Parameters for Acoustic Feedback Cancellation
JP4663630B2 (ja) マルチチャンネルシステム同定装置
Rodríguez et al. Convex Combination of FXECAP–FXECLMS Algorithms for Active Noise Control
Gay Affine projection algorithms
Burra et al. Comparison of convergence performance for LMS and NLMS adaptive algorithms in stereophonic channels
Wahbi et al. Enhancing the quality of voice communications by acoustic noise cancellation (ANC) using a low cost adaptive algorithm based Fast Fourier Transform (FFT) and circular convolution
Schuldt et al. Low-complexity adaptive filtering implementation for acoustic echo cancellation
JP3303898B2 (ja) 適応的伝達関数推定方法及びそれを使った推定装置
Nitsch Real-time implementation of the exact block NLMS algorithm for acoustic echo control in hands-free telephone systems

Legal Events

Date Code Title Description
AK Designated states

Kind code of ref document: A1

Designated state(s): CA JP

AL Designated countries for regional patents

Kind code of ref document: A1

Designated state(s): AT BE CH CY DE DK ES FI FR GB GR IE IT LU MC NL PT SE

ENP Entry into the national phase

Ref document number: 2318929

Country of ref document: CA

WWE Wipo information: entry into national phase

Ref document number: 1999973514

Country of ref document: EP

121 Ep: the epo has been informed by wipo that ep was designated in this application
WWP Wipo information: published in national office

Ref document number: 1999973514

Country of ref document: EP

WWW Wipo information: withdrawn in national office

Ref document number: 1999973514

Country of ref document: EP