Brevet FR2738692A1 — publié le 19970314

PROCÉDÉ DE FILTRAGE NUMÉRIQUE ADAPTATIF DANS LE DOMAINE FRÉQUENTIEL L'invention concerne un procédé de filtrage numérique adaptatif, dans le domaine fréquentiel, de signaux numériques. De tels procédés sont mis en oeuvre en traitement du signal, pour adapter les dispositifs de réception aux évolutions d'un canal de transmission évoluant temporellement. Ainsi, par exemple, en transmission hertzienne, les caractéristiques de propagation des signaux évoluent en fonction des conditions météorologiques: elles sont différentes dans le temps et dans l'espace. On utilise donc généralement, dans les récepteurs, des filtres adaptatifs, c'est à dire des filtres à coefficients de pondération variables dans le temps.Les variations dans le temps des coefficients de pondération sont définies selon un critère d'optimisation et ces coefficients sont produits par des dispositifs (généralement des processeurs de traitement de signal ou DSP) mettant en oeuvre un algorithme d'adaptation. Des coefficients de pondération, calculés par des algorithmes d'adaptation généralement employés, peuvent, dans le domaine temporel, se mettre sous la forme itérative suivante (les indices représentant les dimensions de vecteur ou de matrice) wm(n) = wm(n-l) - 1{J(n)}-G{J(n)} (1) avec n indice temporel, m l'ordre du filtre, y le pas de convergence, J(n) = E(ld(n) - y(n) 12) (fonction de coût minimisée), y(n) les signaux de sortie de filtre correspondants à des signaux d'entrée x(n), d(n) les décisions produites par un organe de décision à partir des signaux de sortie, Hm1{J(n)} une estimation de l'inverse de l'Hessian H de J(n), et est une matrice de dimension m x m. et G{J(n)} une estimation du gradient de J(n) et est un vecteur de dimension m. Cette formulation est très générale et est valable pour une très grande classe d'algorithmes échantillon par échantillon dans le domaine temporel. Les algorithmes correspondant à l'équation (1) sont dits de type quasi-Newton. Un choix connu de G{J(n)} est G{J(n)} = - 2 xm*(n).e(n), avec xm(n) = [x(n),..., x(n - m + 1) jT un vecteur formé des m dernier signaux successifs reçus au temps n, étant le symbole de la conjugaison complexe, T étant une transposition ligne-colonne, et e(n) = d(n) - xmT(n).wm(n - 1) un signal d'erreur. C'est le choix de Hm1{J(n)} qui caractérise principalement les différents algorithmes de quasi Newton. Si cette matrice est la matrice identité, on obtient l'algorithme des moindres carrés moyens (ou LMS, de l'anglais least mean squares). Cet algorithme, très utilisé, présente une complexité faible, au détriment d'une convergence lente également. Si cette matrice est l'inverse de la matrice d'autocorrélation des signaux d'entrée, alors on obtient l'algorithme des moindres carrés récursif (ou RLS, de l'anglais recursive least squares). Cet algorithme présente une convergence meilleure que l'algorithme LMS (les signaux d'entrée étant pris en compte), au prwix d'une complexité beaucoup plus importante. Un problème de mise en oeuvre de cet algorithme d'adaptation est le temps nécessaire à la réalisation du calculs de Hm-l et G pour chaque bloc de signaux d'entrée. D'autres algorithmes de filtrage ont été proposés, dérivés du type RLS, basés sur une extrapolation minmax d'une matrice d'autocorrélation d'un ordre p + 1, inférieur à m, jusqu'à l'ordre m du filtre (G.V. 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). On utilise alors une estimée de l'inverse de la matrice Hessienne donnée par la formule suivante avec j indice variant de O à m - p - 1, avec p Rp-1(n) une matrice carrée de dimension p, p complétée de m-p zéros, b = [0jT:bpH(n - j):l:Om-p-j-1T]H, bmH la transposée hermitienne de bm, bp (t) et ap(t) les prédicteur récursif et puissance d'erreur de prédiction d'ordre p, l'opérateur : représentant une concaténation par juxtaposition, oiT étant un vecteur nul de dimension j transposé, 0m-p-j 1T étant un vecteur transposé d'un vecteur nul de dimension m-p-j-l. Ils permettent d'atteindre une convergence équivalente à un algorithme de type RLS tout en ayant une complexité, proche de celle d'un algorithme de type LMS dans le domaine temporel, nécessitant 2m + 6p multiplications par unité de temps. De tels procédés sont néanmoins de complexité élevée pour des filtres d'ordre élevé, du fait de leur complexité proportionnelle à m (ce qui peut rendre prohibitive leur utilisation pour des applications telles que l'annulation d'échos en audioconférence, où m peut être de l'ordre de plusieurs milliers). Depuis quelques années, dans le domaine du traitement du signal, le filtrage adaptatif dans le domaine fréquentiel a fait l'objet de nombreux développements. Le filtrage adaptatif dans le domaine fréquentiel présente l'avantage, par rapport au filtrage adaptatif dans le domaine temporel, de réduire considérablement la complexité de la mise en oeuvre des algorithmes d'adaptation. Ceci provient du fait que la convolution temporelle est remplacée dans le domaine fréquentiel par une multiplication et deux transformées de Fourier. En conséquence, un calcul dans le domaine fréquentiel permet de réduire à la fois les temps de production des coefficients de pondération et l'encombrement des circuits calculant ces coefficients. Les algorithmes d'adaptation dans le domaine fréquentiel mettent en oeuvre un filtrage adaptatif par blocs, avec une dimension de bloc généralement égale à l'ordre du filtre, c'est à dire au nombre de coefficients de pondération du filtre. La plupart des algorithmes existants sont du type algorithme du gradient. On connaît déjà des exemples de mise en oeuvre d'algorithmes LMS dans le domaine fréquentiel. Ainsi, E.R. Ferrara, "Frequency Domain Adaptive Filtering", Adaptive Filters, C. F. N. Cowan and P. M. Grant, Eds, Englewood Cliffs, Prentice Hall, 1985, Chapitre 6, page 145-179 propose des filtres adaptatifs dans le domaine fréquentiel. Plus particulièrement, la figure 6.2 et la description correspondante présentent un filtre adaptatif de type LMS dans le domaine fréquentiel. En terme de complexité de réalisation, le calcul de m échantillons de sortie d'un filtre à m coefficients avec l'algorithme LMS dans le domaine temporel requiert 2m2 multiplications réelles. Avec l'algorithme décrit dans ce document le même calcul nécessite cinq transformées de Fourier rapides (encore notées FFT) de dimension 2m et deux multiplications complexes sur 2m points. Une FFT sur 2m points demande (m/2) log2m - m multiplications complexes. Le nombre de multiplications complexes par bloc donné est (5m/2) log2m.Une multiplication complexe équivalent à quatre multiplications réelles, le rapport de complexité entre LMS fréquentiel et LMS temporel en terme de multiplications réelles est de (5 log2m + 4)/m. Pour m = 512 par exemple, le rapport sera de 0.096 (ce qui illustre l'intérêt de la mise en oeuvre dans le domaine fréquentiel en terme de complexité de calcul). Un but de l'invention est de proposer un procédé de filtrage adaptatif, dans le domaine fréquentiel, ayant une convergence plus élevée que le filtrage LMS par blocs dans le domaine fréquentiel. Pour ce faire, l'invention propose de modifier le procédé de calcul des blocs de coefficients de pondération en prenant une estimation du Hessian différente de l'identité (c'est à dire à produire à partir des signaux d'entrée un modèle de l'entrée, ce qui permet d'améliorer la convergence), et qui ne nécessite pas le traitement de récursivité inhérent à la formulation de type RLS. Plus précisément, l'invention propose de transposer dans le domaine fréquentiel une version simplifiée du procédé illustré ci-dessus par l'équation (2). L'invention est basée sur la supposition que la matrice d'autocorrélation d'entrée ne change pratiquement pas pendant un intervalle de temps équivalent à au moins m échantillons, c'est à dire qu'une estimée de la matrice d'autocorrélation d'entrée qui est exacte à un temps n sera exacte pour au moins m à 2m intervalles de temps plus tard (un intervalle de temps correspondant à la distance temporelle entre deux signaux d'entrée successifs) . De nombreuses applications répondent à cette hypothèse pour des intervalles de temps beaucoup plus importants (par exemple l'annulation d'échos induit par une propagation en trajets multiples dans la transmission hertzienne). Si on considère que Rm(n) est constante pour au moins m ou 2m intervalles de temps, alors Rp(n) le sera aussi. On peut donc simplifier l'équation (2) de telle sorte que pour j variant de O à (m - 1) - p on a ap(n - j) = ap(n) et bp(n - j) = bp(n). En considérant que G{J(n)} = - 2 xm* (n) e(n), et en écrivant l'équation (1) m fois, avec n = km,..., km + m -1, on obtient une nouvelle équation récursive entre une estimation des coefficients de pondération et l'estimation des coefficients de pondération considérée m intervalles de temps plus tôt wm(k+1) = wm(k) + 24 Zj Rm lxm (j)em j(k) (3), avec j variant de km à km + m - 1, em,j(k) les éléments du bloc d'erreur défini par em(k) = dm(k) - XmT(k)wm(k) (4), et Xm(k) = [xm(km). ..xm(km+m-l)] vecteur d'entrée k étant le rang du bloc. En supposant que la matrice d'autocorrélation d'entrée est invariante pour m intervalles de temps successifs, on en déduit que wm(k + 1) = wm(k) + 2 Rm-1Xrn* (k)em(k) (5), la matrice Rm-1(k) étant celle correspondant au dernier instant (n = km + m - 1) du k-ième bloc. On peut encore écrire Rm-1(k) = Pm(k) + B(k).BH(k) avec Rp -1(k) 0 Pm(k) = 0 0 et Bm(k) est une matrice m * p définie par bp+1 ... m-p-1 Bm(k) = bp+l 0m-p-1 0m-p-2 ... bp+1 avec bp+1 = ap-1/2(k) [bpH(k) 1jH On peut ainsi mettre l'équation (3) sous la forme wm(k + 1) = wm(k) + 2 cm1 (k) + 2 cm2(k) (6) avec cml(k) = Pm(k).Xm*(k).em(k) (7), et cm2(k) = B(k).BH(k).Xm*(k).em(k) (8). On notera que seuls les p premiers éléments de cm1 (k) sont non nuls. Pour le calcul du vecteur cm1 (k) une inversion matricielle et une multiplication matricielle sont nécessaires pour chaque bloc. En choisissant p cml(k) = [(Xp*(k).em(k))H:0m-pH]H avec Xp*(k) une matrice p.m comprenant les p premières rangées de Xm*(k). Autrement dit, le premier terme correctif de l'équation (6) est traité de manière analogue au procédé LMS. Si on fait abstraction du terme correctif cm2 (k) dans l'équation (6), on retrouve une équation de type LMS (à la seule différence que la contrainte de gradient dans le cas de cml(k) est différente du fait que seuls p éléments de ce vecteur sont non nuls). Les convolutions linéaires des équations (4) pour produire le bloc d'erreur et (7) pour le premier bloc correctif peuvent être aisément implantées dans le domaine fréquentiel en utilisant la technique de recouvrement partiel, dite overlap-save, bien connue de l'homme du métier. L'idée de base est de traiter les convolutions linéaires par le biais de convolutions circulaires de taille double, ces convolutions étant mises en oeuvre en utilisant des transformées de Fourier discrète rapide.Ainsi, on obtient c2ml(k) = FQîFlX2mH(k).E2m(k) (9) avec E2m(k) le vecteur bloc d'erreur dans le domaine fréquentiel, F et F-1 l'opérateur de transformée directe de Fourier et son inverse, et Q1 une matrice de dimension 2m x 2m telle que = 1m O Q1 = 00 avec Im la matrice identité de dimension m complété, pour Q1, par des zéros. Si on regarde le second terme correctif de l'équation (6), on peut vérifier d'après l'équation (8) que ce terme a une structure particulière et peut être écrit sous forme d'une succession de trois convolutions linéaires. En terme fréquentiel on pourra se ramener à la formulation suivante C2m2(k) = FQ2F-1B2m*(k)FQ3F-lB2m(k)FQ2F- lX2mH(k)E2m(k) (10) avec B2m(k) définie par B2m(k) = diag{FFT[b*p+1,1:0m-p-1T:b*p+1,p+1 ... b* p+1,2]} avec b* p+1,i le i-ieme élément du vecteur p+l,2 p+l,i b Q2 et Q3 une matrice de dimension 2m x 2m telle que: 1m 0 Im~P 0 Q = et Q = 0 0 0 0 avec 1m-p la matrice identité de dimension m - p. On remarque que l'équation (10) implique cinq transformées de Fourier supplémentaires par rapport à la transposition de l'algorithme LMS dans le domaine fréquentiel, ces transformées étant utilisées pour mettre en oeuvre dans le domaine fréquentiel des contraintes temporelles. Si on considère que p C2m2(k) = B*2m(k) .B2m(k) .X2mH(k).E2m(k) (11). Cette simplification est raisonnable parce que, en supposant que p est très inférieur à m, alors les effets de distorsion sont négligeable et les contraintes matricielles peuvent être supprimées. Si on considère, dans le domaine fréquentiel, l'équation récursive d'adaptation des coefficients de pondération suivante W2m(k+l) = W2m(k) + 2ssC2ml(k) + 2ssC2m2(k) (12) avec C2ml(k) et C2m2(k) tels que définis par les relations (9) et (11), alors on obtient un procédé de filtrage dans le domaine fréquentiel ayant une complexité proche de celle de la transposition fréquentielle de l'algorithme LMS, et présentant une convergence analogue à celle de l'implantation dans le domaine temporel de l'algorithme RLS défini dans la relation (5). Ainsi, l'invention concerne un procédé de filtrage adaptatif de signaux d'entrée temporels, caractérisé en ce que le filtrage est réalisé sur des blocs de m signaux successifs, par multiplication dans le domaine fréquentiel de la transformée de Fourier des blocs de signaux d'entrée par des blocs de coefficients de pondération, ces blocs de coefficients étant calculés récursivement à partir de blocs de coefficients précédents et de premier et deuxième termes correctifs obtenus à partir de signaux de sortie correspondant aux signaux d'entrée filtrés et à partir de la transformée de Fourier d'un modèle de l'entrée. D'autres particularités et avantages apparaîtront à la lecture de la description qui suit, à lire conjointement au dessin annexé dans lequel la figure 1 représente schématiquement un procédé de filtrage adaptatif dans le domaine fréquentiel, selon l'invention. La figure 1 illustre schématiquement un exemple de mise en oeuvre du procédé de filtrage adaptatif selon l'invention. L'objet du procédé est de pondérer des signaux d'entrée x(n) successifs par des coefficients wi(k) variables et de produire des signaux de sortie y(n) correspondants tels qu'on minimise des signaux d'erreur e(n) = d(n) - y(n), avec d(n) des décisions produites par un organe de décision à partir des signaux de sortie y(n). Une application d'un tel procédé est par exemple l'annulation d'échos en audioconférence. Une autre application possible est l'identification adaptative d'une manière générale, et en particulier l'identification de canaux perturbés par des échos lors de transmissions ou de diffusions hertziennes. Dans ce contexte, l'ordre des prédicteurs ainsi que les conditions de stationnarité rentrent exactement dans les hypothèses de l'invention. Dans l'invention, les signaux d'entrée x(n) sont traités par blocs de m signaux successifs (x(n),.... x(n+m-1)), avec m entier. Généralement m est une puissance de 2 (par exemple m = 512 ou 1024). Classiquement, à un bloc de signaux d'entrée (x(n) ,.. . , x(n+m-l)) va correspondre un bloc de m coefficients de pondération (wn(k),..., Wn+m-l(k)). Schématiquement, le procédé selon l'invention propose de réaliser, dans le domaine fréquentiel, l'adaptation des coefficients d'une part à partir d'un bloc d'erreurs produit à partir des signaux de sortie (de manière analogue à la réalisation dans le domaine fréquentiel d'un procédé de type LMS) et d'autre part à partir d'un prédicteur de signaux d'entrée. On va décrire maintenant plus précisément le procédé illustré sur la figure. Les signaux d'entrée x(n) sont regroupés, par une transformation série-parallèle 1, en blocs de signaux (x(km),..., x(km+m-l)). Chaque k-ième bloc xm(k) est regroupé par une concaténation 2, avec le bloc précédent xm(k-l) , puis le bloc de 2m signaux obtenus est transposé dans le domaine fréquentiel, par une transformée de Fourier rapide 3 (notée FFT). On obtient alors un bloc fréquentiel de 2m signaux d'entrée. Par une multiplication 4, on produit 2m signaux représentatifs du résultat de l'opération suivante Y2m(k) = X2m(k)W2m(k) avec X2m(k) = diag(FFT[x(km-m)....x(km+m-l)lT) et W2m(k) un vecteur, de dimension 2m, de coefficients de pondération. Après une transformation de Fourier rapide inverse 5 de Y2m(k) et une déconcaténation 6 on produit un bloc y, (k) de m signaux correspondant à la seconde moitié des signaux obtenus après la transformation inverse, dans le domaine temporel. Autrement dit, le bloc y,(k) est représentatif de la relation suivante y,(k) = Q0F-ly2m(k) avec = = [0m, Im] une matrice de dimension m.2m, avec 0m une matrice carrée nulle de dimension m, 1m la matrice identité de dimension m, et F-1 l'inverse de l'opérateur de transformée directe de Fourier de dimension 2m. Par une transformation parallèle-série 7, on obtient les signaux de sortie y(n). En ce qui concerne l'adaptation des coefficients de pondération, on réalise l'opération suivante W2m(k+l) = W2m(k) + 2bC2ml(k) + 2ssC2m2(k) (12), avec C2m1(k) = FQ1F-1x2mH(k)#E2m(k) (9) et C2m2(k) = B2m*(k).B2m(k).X2mH(k).E2m(k) (11) avec un pas de convergence, E2m(k) un vecteur bloc d'erreur dans le domaine fréquentiel, de dimension 2m, X2mH(k) une matrice diagonale de dimension 2m.2m obtenue par une transposition hermitienne 8 de X2m(k), F l'opérateur de transformée directe de Fourier de dimension 2m, Q1 une matrice de dimension 2m * 2m telle que p entier inférieur à m, Ip matrice identité de dimension p, B2m(k) matrice de dimension 2m définie par B2m(k) = diag{FFT[b*p+1,1:Omp1T:: b*p+î,p+î b*p+1,2]T} avec b*p+1 i le i-ième élément d'un vecteur b*p+1 de dimension p+l. Pour produire c2m1(k), on procède de la manière suivante - à partir du bloc ym(k) de m signaux de sortie, définis dans le domaine temporel, et d'un bloc dm(k) de m décisions correspondantes, on produit par une soustraction 9 un bloc de m signaux d'erreur em(k) = dm(k) - ym(k). Les décisions sont produites par un circuit de décision non linéaire, - après une concaténation 10 de em(k) avec un bloc de m signaux nuls, on produit le bloc E2m(k) de 2m signaux d'erreurs dans le domaine fréquentiel par une transformée de Fourier rapide 11 (on a donc E2m(k) = FQOTem(k)), - à partir de E2m(k) on produit, par une multiplication 12 avec le bloc X2mH(k) de 2m signaux, un bloc de 2m signaux représentatifs de l'opération X2mH(k) E2m(k), - ce bloc de 2m signaux obtenu est ensuite transposé dans le domaine temporel par transformée de Fourier inverse rapide 13, les m derniers signaux correspondants sont remplacés par des zéros lors d'une concaténation 14, et les 2m signaux ainsi produits sont transposés dans le domaine fréquentiel par une transformée de Fourier rapide 15. Autrement dit, le terme correctif C2m1(k) est obtenu par un procédé semblable à celui utilisé dans la mise en oeuvre d'un algorithme de type LMS dans le domaine fréquentiel, en remplaçant 12m par Q1. Le terme correctif C2m2(k) correspond à une modélisation des signaux d'entrée. L'invention est basée sur la supposition que la matrice d'autocorrélation d'entrée ne change pratiquement pas pendant un intervalle de temps équivalent à au moins m échantillons, c'est à dire qu'une estimée de la matrice d'autocorrélation d'entrée qui est exacte à un temps n sera exacte pour au moins m à 2m intervalles de temps plus tard (un intervalle de temps correspondant à la durée entre deux signaux d'entrée successifs), voire plus. De nombreuses applications répondent à cette hypothèse pour nombre l.m d'intervalles de temps beaucoup plus important (par exemple l'annulation d'échos induits par une propagation en trajets multiples dans la transmission hertzienne). En conséquence, on considère des regroupements dans un concaténateur 16 de 1 blocs successifs de signaux d'entrée. Pour chaque regroupement, on calcule dans le domaine temporel un bloc bp+1 de p+l signaux (modèle de l'entrée) tels que bp+1 = ap-1/2(k) [-bpH(k) 1]H avec bp et ap des prédicteur d'entrée et puissance d'erreur de prédiction d'ordre p, avec p Pour ce faire, on utilise dans un corrélateur 17 la méthode connue dite d'autocorrélation (décrite par exemple dans "Digital Spectral Analysis with Applications", par S.L. Marple, Prentice Hall, New Jersey, 1987). A partir du bloc bp+1 on produit par une concaténation 18 dans le domaine temporel le bloc b'm de m signaux (b*p+1,1 0m-p-1T b*+1,p+1 ... b* p+1,2). On passe ensuite dans le domaine fréquentiel par une transformée de Fourier rapide 19 pour obtenir un bloc de 2m échantillons représentatif de B2m(k) tel que: B2m (k)-diag{FFT[b*p+1,1:0m-p-I T:b* p+1,o+1...b*p+1,2]T}. Pour obtenir C2m2(k), on procède ensuite à la multiplication 20 de B2m(k) par sa conjuguée B2m*(k) (obtenue par une transformation 21), et à la multiplication 22 par X2mH(k).E2m(k), ce terme ayant été calculé par ailleurs pour produire c2ml(k). A partir des deux blocs c2ml(k) et C2m2 (k) de 2m signaux correctifs, et après multiplications 22 et 23 de ces blocs par 2ss, on peut calculer pour chaque bloc de signaux d'entrée le bloc de coefficients de pondération correspondant par une simple addition 24 des signaux du bloc de coefficients de pondération précédent fourni par un circuit de retard 25 avec les signaux des blocs 2ssC2m1(k) et 2ssC2m2(k). REVENDICATIONS 1 - Procédé de filtrage adaptatif de signaux d'entrée (x(n)) temporels, caractérisé en ce que le filtrage est réalisé sur des blocs (xm(k)) de m signaux successifs ((x(n),...,x(n+m-1))), par multiplication (4) dans le domaine fréquentiel de la transformée de Fourier (X2m(k)) des blocs de signaux d'entrée par des blocs (W2m(k)) de coefficients de pondération, ces blocs de coefficients étant calculés récursivement (915) à partir de blocs de coefficients précédents et de premier et deuxième termes correctifs (C2m1(k), C2m2(k)) obtenus d'une part (13-15) à partir de signaux de sortie (ym(k)) correspondants aux signaux d'entrée filtrés et d'autre part (16-20) à partir de la transformée de Fourier d'un modèle de l'entrée (bp+1). 2 - Procédé selon la revendication 1, caractérisé en ce que le modèle de l'entrée (bop+1) est calculé par autocorrélation (17) et comprend un nombre p de signaux, avec p inférieur à m. 3 - Procédé selon la revendication 2, caractérisé en ce qu'un modèle de l'entrée est calculé (16) pour l.m blocs de signaux d'entrée successifs. 4 - Procédé selon l'une des revendications 1 à 3, caractérisée en ce que le calcul du second terme correctif (C2m2(k)) comprend la transformation par transformée de Fourier (19) d'un bloc (b'm) de m signaux produits dans le domaine temporel à partir du modèle de l'entrée (bop+1). 5 - Procédé selon l'une des revendications 1 à 4, caractérisé en ce que le premier terme correctif (C2ml(k)) est obtenu par la mise en oeuvre (10-15) dans le domaine fréquentiel de l'algorithme des moindres carrés moyens.