EP0848867A1

EP0848867A1 - Adaptive digital filtering process in the frequency domain

Info

Publication number: EP0848867A1
Application number: EP96931093A
Authority: EP
Inventors: Constantinos Berberidis; Jacques Palicot
Original assignee: Telediffusion de France ets Public de Diffusion; France Telecom SA
Current assignee: Telediffusion de France ets Public de Diffusion; Orange SA
Priority date: 1995-09-08
Filing date: 1996-09-09
Publication date: 1998-06-24
Also published as: FR2738692B1; FR2738692A1; WO1997009782A1

Abstract

The invention concerns an adaptive filtering process for time input signals (X(n)) in the frequency domain. Filtering is carried out on blocks (xm(k)) of m successive signals ((x(n)), ... x (n+m-1))), by multiplying the Fourier transform (X2m(k)) of the input signal blocks by blocks (W2m(k)) of weighting coefficients in the frequency domain, the coefficient blocks being calculated recursively from preceding coefficient blocks and from first and second corrective terms (C2m1(k), C¿2m?2(k)) obtained from output signals (y¿m?(k)) corresponding to the filtered input signals and from the Fourier transform of a model of the input (bp+1).

Description

Adaptive digital filtering process in the frequency domain.

The invention relates to a method of adaptive digital filtering, in the frequency domain, of digital signals. Such methods are implemented in signal processing, to adapt the reception devices to changes in a time-changing transmission channel. Thus, for example, in hertzian transmission, the propagation characteristics of signals evolve according to weather conditions: they are different in time and in space. Adaptive filters are therefore generally used in receivers, that is to say filters with weighting coefficients that vary over time. The variations in time of the weighting coefficients are defined according to an optimization criterion and these coefficients are produced by devices (generally signal processing processors or DSPs) implementing an adaptation algorithm.

Weighting coefficients, calculated by generally used adaptation algorithms, can, in the time domain, take the following iterative form (the indices representing the vector or matrix dimensions): ^w m < ⁿ ) ^{= w} m ( ^{n_: L} ) "^ _m ^_1 {J (n)} G {J (n) ^} (1) with: n time index, m order of the filter, μ the convergence step,

J (n) = E {| d (n) - y (n) | ² } (minimized cost function;) y (n) the corresponding filter output signals to input signals x (n), d (n) the decisions produced by a decision-making body from the output signals,

H_ ^-1 {J (n)} an estimate of the Hessian inverse (H) of J (n) and is a matrix of dimension mxm, and G {J (n)} an estimate of the gradient of J ( n) and is a vector of dimension m.

This formulation is very general and is valid for a very large class of sample by sample algorithms in the time domain. The algorithms corresponding to equation (1) are said to be of the quasi-Newton type.

A known choice of G {J (n)} is:

G {J (n)} = - 2 x _m * (n) e (n), with ^x m ( ⁿ ) ⁼ [ ^χ ( ⁿ ) '- * -' ^χ ( ⁿ - m + 1)] ^τ a vector formed of the last m successive signals received at time n, * being the symbol of the complex conjugation, T being a row-column transposition, and e (n) = d (n) - x _m ^τ (n) w _m (n - 1) an error signal. It is the choice of # _m ^_1 {J (n)} which mainly characterizes the different quasi- Newton algorithms.

If this matrix is the identity matrix, we obtain the least mean squares algorithm (or LMS). This algorithm, widely used, has low complexity, at the expense of slow convergence as well. For example, the document IEEE ACOUSTICS, SPEECH, AND SIGNAL PROCESSING MAGAZINE Vol. 34 No., February 1986; New York, US, pages 105-117 entitled "Block realization of multirate adaptive digital filters" proposes a block formulation of the LMS algorithm in the time domain.

If this matrix is the inverse of the autocorrelation matrix of the input signals, then we obtains the recursive least squares algorithm (or RLS, from English recursive least squares). This algorithm presents a better convergence than the LMS algorithm (the input signals being taken into account), at the cost of a much greater complexity. A problem with the implementation of this adaptation algorithm is the time necessary for carrying out the calculations of tf _m ^-1 and G for each block of input signals. An example of implementation of this algorithm is described in document US-A-4 658 426. Another document IEEE, Acoustics Speech and Signal Processing Magazine, Vol. 34 No. 6 December 1986, New York, US, pages 1573-1582, entitled "Self orthogonalizing efficient block adaptive filters" proposes a block-type RLS algorithm in the time domain.

Other filtering algorithms have been proposed, derived from the RLS type, based on a min-max extrapolation of an autocorrelation matrix of order p + 1, less than m, up to order m of the filter. (GV Moustakides, S. Theodoridis, "Fast Quasi Newton Transversal Filters - A New Class of Adaptive Estimation Algorithms", IEEE Trans. On Signal Processing, vol. 39, pp 2184-2193, Oct. 1991). We then use an estimate of the inverse of the Hessian matrix given by the following formula:

* V ^> m ^H (2) with j index varying from O to m - p - 1, with p <m, R_ ₎ ~ ¹ (n) a square matrix of dimension p, completed with mp zeros, b _m = _[0j T _: - _bp H (n - j): 1: O ^^^ r-, bm ^H la h-xev-ru.m ... i ... tV.i._eWn „nΛ.e ~ d ~ e_ b ~ _ _m ,, b (t) and a (t) the recursive predictor and power ir xr of order prediction error p. the operator: representing a concatenation by juxtaposition, 0- ^τ being a zero vector of dimension j transposed, O _m _ _; _ ₁ ^τ being a vector transposed from a zero vector of dimension mpj-1.

They allow to reach a convergence equivalent to an RLS type algorithm while having a complexity, close to that of an LMS type algorithm in the time domain, requiring 2m + 6p multiplications per unit of time. Such methods are nevertheless of high complexity for filters of high order, because of their complexity proportional to m (which can make their use prohibitive for applications such as the cancellation of echoes in audio conference, where m can be several thousand).

In recent years, in the field of signal processing, adaptive filtering in the frequency domain has been the subject of numerous developments. Adaptive filtering in the frequency domain has the advantage, compared to adaptive filtering in the time domain, of considerably reducing the complexity of the implementation of adaptation algorithms. This comes from the fact that the temporal convolution is replaced in the frequency domain by a multiplication and two Fourier transforms. Consequently, a calculation in the frequency domain makes it possible to reduce both the production times of the weighting coefficients and the size of the circuits calculating these coefficients.

The frequency domain adaptation algorithms implement adaptive block filtering, with a block dimension generally equal to the order of the filter, i.e. the number of filter weights. Most of the existing algorithms are of the gradient algorithm type. Examples of implementation of LMS algorithms in the frequency domain are already known. ER Ferrara, "Frequency Domain Adaptive Filtering", Adaptive Filters, CFN Cowan and PM Grant, Eds, Englewood Cliffs, Prentice Hall, 1985, Chapter 6, pages 145-179 offers adaptive filters in the frequency domain. More particularly, figure 6.2 and the corresponding description present an adaptive filter of the LMS type in the frequency domain. In terms of complexity of realization, the calculation of m output samples of a filter with m coefficients with the LMS algorithm in the time domain requires 2m ² real multiplications. With the algorithm described in this document the same calculation requires five fast Fourier transforms (also noted FFT) of dimension 2m and two complex multiplications on 2m points. An FFT 2m application point (m / 2) log ₂ m - m complex multiplications. The number of complex multiplications is given block (5m / 2) log ₂ m. A complex multiplication equivalent to four real multiplications, the complexity ratio between frequency LMS and temporal LMS in terms of real multiplications is (5 log ₂ m + 4) / m. For m = 512 for example, the ratio will be 0.096 (which illustrates the interest of the implementation in the frequency domain in terms of computational complexity). An object of the invention is to propose an adaptive filtering method, in the frequency domain, but having a higher convergence than LMS block filtering in the frequency domain. To do this, the invention proposes to modify the calculation method blocks of weighting coefficients by taking a Hessian estimate different from the identity (that is to say from producing the input signals a model of the input, which makes it possible to improve convergence), and which does not require the recursion processing inherent in the RLS type formulation.

More specifically, the invention proposes to transpose into the frequency domain a simplified version of the method illustrated above by equation (2). The invention is based on the assumption that the input autocorrelation matrix practically does not change during a time interval equivalent to at least m samples, that is to say that an estimate of the autocorrelation matrix of input which is exact at a time n will be exact for at least m to 2 m time intervals later (a time interval corresponding to the time distance between two successive input signals). Many applications meet this hypothesis for much longer time intervals (for example the cancellation of echoes induced by a multipath propagation in the radio transmission).

If we consider that R ^ n) is constant for at least m or 2m time intervals, then R_ (n) will be also.

We can therefore simplify equation (2) so that for j varying from 0 to (m - 1) - p we have: a _p (n - j) = a _p (n) and b _p (n - j ) = b _p (n). Considering that G {J (n)} = - 2 x _m * (n). e (n), and by writing equation (1) m times, with n = km, ..., km + m -1, we obtain a new recursive equation between an estimate of the weights and the estimate of the weighting coefficients considered m earlier time intervals: w _m (k + l) = w _m (k) + 2μ ∑j R _m - ¹ x _m * (j) e _{m ι} j (k) (3), with j varying from km to km + m - 1, e _m -: (k) the elements of the error block defined by: e _m (k) = d _m (k) - X _m T (k) w _m (k) (4), and ^x m ^k ^ ⁼ [ ^χ m ( ^km ) * • • • ^x _m ( ^{kln + m_1} ) 1 input vector k being the rank of the block.

Assuming that the input autocorrelation matrix is invariant for m successive time intervals, we deduce that: w _χn (k + 1) = w _m (k) + 2β R _m - ¹ X _m * (k) e _ιn (k) (5), the matrix R ^ ^-1 (k) being that corresponding to the last instant (n = km + m - 1) of the k-th block. We can still write:

and B _m (k) is a matrix m * p defined by

We can thus put equation (3) in the form: w _m (k + 1) = w _m (k) + 2μ ^ (k) + 2μ c _m ² (k) (6) with c _m l (k ) = P _m (k) X _m * (k) _eιn (k) (7), and c _m ² (k) = B (k) BH (k) X _m * (k) _eιn (k) (8) .

Note that only the first p elements of c _m ¹ (k) are non-zero. For the calculation of the vector c _m ¹ (k) a matrix inversion and a matrix multiplication are necessary for each block. By choosing p << m we can choose: c _m l (k) = [(X _p * (k) .e _m (k)) H: O _m _ _p H] H with X _p * (k) a matrix p. m including the first p rows of X _m * (k).

In other words, the first corrective term of equation (6) is treated in a similar way to the LMS process. If we ignore the corrective term c _m ² (k) in equation (6), we find an equation of LMS type (with the only difference that the gradient constraint in the case of c _m ¹ (k) is different because only p elements of this vector are non-zero).

The linear convolutions of equations (4) to produce the error block and (7) for the first corrective block can be easily implemented in the frequency domain using the technique of recovery of partial recovery called overlap-save, well known in the art. skilled in the art. The basic idea is to treat linear convolutions by means of circular convolutions of double size, these convolutions being implemented using fast discrete Fourier transforms. Thus, we obtain: ^ m ¹ **) ^{≈ F} Q ¹ F- ¹ X _2m ^H (k) .E _2ιn (k) (9) with

E _2m (k) the error block vector in the frequency domain,

F and F " ¹ the direct Fourier transform operator and its inverse, and Q ¹ a matrix of dimension 2m x 2m such that:

with l _m the identity matrix of dimension m, completed for Q ¹ by zeros.

If we look at the second corrective term of equation (6), we can verify from equation (8) that this term has a particular structure and can be written in the form of a succession of three linear convolutions. In frequency terms, we can bring back the following formulation:

C _2m ² (k) = FQ2 _F -lB _2m * (k) FQ3 _F - ^l B _2ιn (k) FQ ² _F -lχ _2m H (k) E _2m (k) (10) with

B _2m (k) defined by:

B _2m (k) = diag {FFT [b * _{p + 1/1} 0 _m _ _p _ ^ b * p + l, p + l "* Vl, ^{2] T} > with p + l, i the i-th element of the vector b * p + 1

Q ² and O ³ a matrix of dimension 2m x 2m such that:

Q ² = ^L m and _Q 3 = ^L m

0 0 with I -, _m _ _D the identity matrix of dimension m - p.

Note that equation (10) implies five additional Fourier transforms compared to the transposition of the LMS algorithm in the frequency domain, these transforms being used to implement temporal constraints in the frequency domain.

If we consider that p << m, then we can simplify the expression of C ₂ ² (k) so that:

C ° 2m ² (k) = B ¹ 2m (k) B 2m '(k) X _2m ^H (k) E 2m (k) (11). This simplification is reasonable because, assuming that p is much less than m, then the distortion effects are negligible and the matrix constraints can be eliminated.

If we consider, in the frequency domain, the recursive equation of adaptation of the following weighting coefficients:

W 2m '(k + 1) = W _{? M} (k) + 2μC' 2ml (k) + 2μC 2m ² (k) (12) with

C _{? M} ¹ (k) and C '2-> m ² (k) as defined by the relations

(9) and (11), then a filtering process in the frequency domain is obtained having a complexity close to that of the frequency transposition of the LMS algorithm, and having a convergence analogous to that of the implementation in the time domain of the RLS algorithm defined in relation (5).

Thus, the invention relates to an adaptive filtering method of time input signals, in which the filtering is carried out on blocks of m successive signals, by multiplication in the frequency domain of the Fourier transform of blocks of input signals by blocks of weighting coefficients, these blocks of coefficients being computed recursively from blocks of previous coefficients and from first and second corrective terms obtained from output signals corresponding to the filtered input signals and from the transform of Fourier of an input model, characterized in that the input model is calculated by autocorrelation.

Other particularities and advantages will appear on reading the following description, to be read in conjunction with the appended drawing in which FIG. 1 schematically represents an adaptive filtering process in the frequency domain, according to the invention.

FIG. 1 schematically illustrates an example of implementation of the adaptive filtering method according to the invention.

The object of the method is to weight successive input signals x (n) by variable coefficients wWk) and to produce corresponding output signals y (n) such that error signals e (n) are minimized = d (n) - y (n), with d (n) decisions produced by a decision-making body from the output signals y (n).

An application of such a method is for example the cancellation of echoes in audioconference. Another possible application is adaptive identification in general, and in particular the identification of channels disturbed by echoes during radio transmissions. In this context, the order of the predictors as well as the stationarity conditions fall exactly within the hypotheses of the invention.

In the invention, the input signals x (n) are processed in blocks of m successive signals (x (n), ..., x (n + m-l)), with m integer. Generally m is a power of 2 (for example m = 512 or 1024).

Conventionally, a block k of input signals (x (n), ..., x (n + ml)) will correspond to a block of m weighting coefficients (w _n (k), ..., w _{n + m} _ ₁ (k)).

Schematically, the method according to the invention proposes to carry out, in the frequency domain, the adaptation of the coefficients on the one hand from a block of errors produced from the output signals (analogously to the realization in the frequency domain of an LMS type method) and on the other hand from a predictor of input signals.

We will now describe more precisely the process illustrated in the figure. The input signals x (n) are grouped, by a series-parallel transformation 1, into signal blocks [x (km), ..., x (km + m-l)].

Each k-th block ^x _m (k) ^is grouped by a concatenation 2, with the previous block x _ιn (kl), then the block of 2m signals obtained is transposed in the frequency domain, by a fast Fourier transform 3 (denoted FFT). We then obtain a frequency block of 2m input signals.

By a multiplication 4, we produce 2m signals representative of the result of the following operation:

^Y 2m ( ^k > = ^X 2m ( ^k > ^W 2m ( ^k ) ^with

X _2m (k) = diag {FFT [x (km-m) .... x (km + ml)] ^τ } and

W _2m (k) a vector, of dimension 2m, of weighting coefficients.

After a fast inverse Fourier transformation of Y2m ( ^k ) ^e * " ^a deconcantenation 6 we produce a block y _m (k) of m signals corresponding to the second half of the signals obtained after the inverse transformation, in the time domain.

In other words, the block y _m (k) is representative of the following relation: y _m (k) = Q ° F- ¹ Y _2m (k) with

Q ° = [0 _m , I _m ] a matrix of dimension m. 2m, with 0 _m a zero square matrix of dimension m, I the identity matrix of dimension m, and F ^-1 the inverse of

The 2m direct Fourier transform operator.

By a parallel-series transformation 7, the output signals y (n) are obtained.

As regards the adaptation of the weighting coefficients, the following operation is carried out:

W _2m (k + 1) = W _2m (k) + 2μC _2m ¹ (k) + 2μC _2m ² (k) (12), with

^ m ¹ ^) ^{= FQlF} ~ ^lχ 2m ^H ( ^k ) ' ^E 2m ( ^k) < ⁹ > ^and C _2m ² (k) = B _2m * (k) .B _2m (k) .X _2ln H _(k) .E _2ln (k) (11) with μ a convergence step,

E _2m (k) an error block vector in the frequency domain, of dimension 2m, ^X 2m ^H ( ^k ) ^a diagonal matrix of dimension 2m.2m obtained by Hermitian transposition 8 of X2m ( ^k ) '

F the direct Fourier transform operator of dimension 2m, Q ¹ a matrix of dimension 2m * 2m such that:

* P °

0 0 p integer less than m, I _D identity matrix of dimension p, B _2m (k) matrix of dimension 2m defined by:

- from the block y _ιn (k) of m output signals, defined in the time domain, and from a block d _m (k) of m corresponding decisions, a block of m signals of d is produced by subtraction 9 error e _m (k) = d _m (k) - y _ιrι (k). Decisions are produced by a non-linear decision circuit,

- after a concatenation 10 of e _m (k) with a block of m null signals, we produce the block E _2m (k) of 2m error signals in the frequency domain by a fast Fourier transform 11 (we therefore have E _2ιn (k) = FQθT _eιn (k)), from E _2m (k) we produce, by a multiplication 12 with the block X _2m ^H (k) of 2m signals, a block of 2m signals representative of the operation

^X 2m ^H ( ^k > - ^E 2m ( ^k ) <

- this block of 2m signals obtained is then transposed in the time domain by fast inverse Fourier transform 13, the last m corresponding signals are replaced by zeros during a concatenation 14, and the 2m signals thus produced are transposed in the domain frequency by a fast Fourier transform 15. In other words, the corrective term C _2m ¹ (k) is obtained by a process similar to that used in the implementation of an LMS type algorithm in the frequency domain, by replacing I _2m by Q ¹ .

The corrective term ^C ₂ m ² ( ^k ) corresponds to a modeling of the input signals.

The invention is based on the assumption that the input autocorrelation matrix practically does not change during a time interval equivalent to at least m samples, that is to say that an estimate of the autocorrelation matrix of entry which is exact at a time n will be exact for at least m to 2. m later time intervals (a time interval corresponding to the duration between two successive input signals), or even more. Many applications meet this hypothesis for number 1. m much longer time intervals (for example the cancellation of echoes induced by a multipath propagation in the hertzian transmission). Consequently, we consider groupings in a concatenator 16 of 1 successive blocks of input signals. For each grouping, we calculate in the time domain a block b _{p + 1} of p + 1 signals (input model) such that b _{p + 1} = a _p -l / ² (k) [-b _p H ( k) 1] ^H with b _D and a _n of the input predictors and prediction error power of order p, with p <m.

To do this, a known method known as autocorrelation is used in a correlator 17 (described for example in "Digital Spectral Analysis with Applications", by S.L. Marple, Prentice Hall, New Jersey, 1987).

From block b> ₊₁ - block b 'is produced by a concatenation 18 in the time domain of m signals (b * _{+ 1 # 1} , 0 ^. ^, b * _{+ lf p + 1} , ..., b * _{p + 1 ^ 2} ).

We then go into the frequency domain by a fast Fourier transform 19 to obtain a block of 2m samples representative of B _2m (k) such that: B 2, _m m (* k) '= diag ^ { ^l FFT [ ^l b ^* p + .1.,. 1: 0 _m m- „p-1. ^τ : b ^* p + l - ,, p + .l .... b ^* p + .l., 2 ₀ ] ^JT } '.

To obtain ^c 2m ² ( ^k ) ' ^we then proceed to the multiplication 20 of B _2m (k) by its conjugate B ₂ _ * (k) (obtained by a transformation 21) and to the multiplication 22 by X _2m ^H (k ) E _2m (k), this term having been calculated elsewhere to produce C _2m ¹ (k).

From the two blocks C _2m ¹ (k) and C _2m ² (k) of 2m corrective signals, and after multiplication 22 and 23 of these blocks by 2μ, one can calculate for each block of input signals the block of coefficients of corresponding weighting by a simple addition 24 of the signals of the block of preceding weighting coefficients, supplied by the delay circuit 25, with the signals of the blocks 2μC _2m l (k) and 2μC _2m ² _(k) .

Claims

1 - Method for adaptive filtering of time input signals (x (n)), in which the filtering is carried out on blocks (x _m (k)) of m successive signals ((x (n), ..., x (n + ml))), by multiplication (4) in the frequency domain of the Fourier transform (3) (X _2ln (k)) of the input signal blocks by blocks (W _2ιn (k)) of weighting coefficients, these blocks of coefficients being calculated recursively (9-15) from blocks of previous coefficients and of first and second corrective terms ( ^C 2m ^{1 k} ) 'C _2jn ² (] _C )) obtained from a part (13-15) from output signals (y _rn (k)) corresponding to the filtered input signals and secondly (16-20) from the Fourier transform of a model of the input (b _{p + 1} ), characterized in that the input model (b _{p + 1} ) is calculated by autocorrelation (17) of the input signals.

2 - Method according to claim 1 characterized in that it comprises a number p of signals, with p less than m. 3 - Method according to claim 2, characterized in that a model of the input is calculated for (16) 1. m blocks of successive input signals.

4 - Method according to one of claims 1 to 3, characterized in that the calculation of the second corrective term ( ^c 2m ² ^ ^k ^ comprises the transformation by Fourier transform (19) of a block ( ^b, _m ) ^of m signals produced in the time domain from the input model (b _{p + 1} ).

5 - Method according to one of claims 1 to 4, characterized in that the first corrective term

(c _2m ¹ (k)) is obtained by the implementation in the frequency domain of 1 * -a: the least squares algorithm means .