Brevet FR2888963A1 — publié le 20070126

La présente invention concerne les procédés de poursuite d'un espace de dimension r inférieure à la dimension n des vecteurs de données, notamment dans le cadre de l'analyse spectrale de données audio. L'invention concerne plus particulièrement les procédés dits d'analyse haute résolution (HR). Le calcul de l'espace signal peut permettre par exemple de décomposer des données observées en une somme d'exponentielles complexes. La dimension r est par exemple le nombre d'exponentielles que l'on choisit pour modéliser le signal utile. L'analyse d'un signal en termes d'espaces est une technique désormais classique de traitement du signal, dans laquelle l'espace signal contenant toute l'information caractérisant le signal est distingué de l'espace bruit qui ne contient que du bruit. Cette approche a été principalement étudiée dans les domaines du filtrage adaptatif, de la localisation de sources et de l'estimation paramétrique, notamment dans la publication [P. COMON et G.H. GOLUB, Tracking a few extreme singular values and vectors in signal processing , Proc. of IEEE, volume 78, pages 1327-1343, août 1990]. La décomposition en valeurs propres (EVD) et la décomposition en valeurs singulières (SVD) sont des outils classiques pour estimer l'espace signal. Leur principal inconvénient est leur complexité algorithmique élevée. C'est pourquoi, dans le cadre d'un traitement adaptatif du signal, il existe un intérêt réel à développer des techniques rapides de poursuite de l'espace signal. Parmi les techniques rapides proposées à ce jour, figure l'algorithme décrit dans l'article [I. KARASALO, Estimating the covariance matrix by signal subspace averaging , IEEE Trans. Acoust., Speech, Signal Processing, 34:8-12, février 1986], qui approche la décomposition en valeurs singulières de la matrice de données en calculant la décomposition en valeurs singulières d'une matrice plus petite. Une approche similaire a été développée dans [E. C. REAL, D. W. TUFTS et J.W. COOLEY, Two algorithms for fast approximate subspace tracking , IEEE Trans. Signal Processing, 47(7) :1936-1945, juillet 1999]. L'algorithme FST présenté dans [D.J. RABIDEAU, Fast, rank adaptive subspace tracking and applications , IEEE Trans. Signal Processing, 44(9) :2229-2244, septembre 1996] remplace la décomposition en valeurs singulières décrite dans [I. KARASALO, Estimating the covariance matrix by signal subspace averaging , IEEE Trans. Acoust., Speech, Signal Processing, 34:8-12, février 1986] par des rotations de Givens, conduisant à un suivi plus rapide. Une autre approche consiste à entrelacer une mise à jour récursive de la matrice de corrélation ou de la matrice de données avec une ou plusieurs étapes d'un algorithme standard de décomposition en valeurs singulières. C'est le cas de la méthode de Jacobi décrite dans [M. MOONEN, P. Van DOOREN J. VANDEWALLE, An SVD updating algorithm for subspace tracking , SIAMJ. Matrix Ana. Appl., 13(4):1015-1038, 1992], de l'itération QR transposée décrite dans [E. M. DOWLING, L. P. AMMANN R. D. DEGROAT, A TQR-iteration based adaptive SVD for real time angle and frequency tracking , IEEE Trans. Signal Processing, 42(4):914-926, 1994], de l'itération orthogonale/bi-orthogonale décrite dans [R. BADEAU, G. RICHARD B. DAVID, Sliding window adaptive SVD algorithms , IEEE Trans. Signal Processing, 52(1):1-10, 2004], dans [P. STROBACH, Low-rank adaptive filters , IEEE Transactions on Signal Processing, 44(12):2932- 2947, 1996], dans [P. STROBACH, Bi-iteration SVD subspace tracking algorithms , IEEE Trans. Signal Processing, 45(5):1222-1240, 1997], et dans [P. STROBACH, Square hankel SVD subspace tracking algorithms , Signal Processing, 57(1):1-18, 1997], et de la méthode des puissances itérées décrite dans [Y. HUA, Y. XIANG, T. CHEN, K. ABEDMERAIM Y. MIA0, A new look at the power method for fast subspace tracking , Digital Signal Processing, 1999]. Certaines techniques de poursuite sont fondées sur d'autres décompositions matricielles, comme la factorisation QR révélatrice de rang décrite dans [C. H. BISCHOF G. M. SHROFF, On updating signal subspaces , IEEE Trans. Signal Processing, 40:96-105, 1992], la décomposition URV révélatrice de rang décrite dans [G. W. STEWART, An updating algorithm for subspace tracking , IEEE Trans. Signal Processing, 40:1535-1541, 1992], et la (bi)-diagonalisation de Lankzos décrite dans [G. XU, H. ZHA, G. H. GOLUB T. KAILATH, Fast algorithms for updating signal subspaces , IEEE Trans. Circuits Syst., 41(8):537-549, 1994]. Une approche conceptuellement différente consistant à considérer l'estimation de l'espace signal comme un problème d'optimisation avec ou sans contraintes, est retenue dans [T. CHEN S. AMARI, Unified stabilization approach to principal and minor components extraction algorithms , Neural Networks, 14(10):1377-1387, 2001], dans [Z. Fu E. M. DOWLING, Conjugate gradient eigenstructure tracking for adaptive spectral estimation , IEEE Trans. Signal Processing, 43(5):1151-1160, 1995], dans [S. Y. KUNG, K. I. DIAMANTARAS J. S. TAUR, Adaptive principal component extraction (apex) and applications , IEEE Trans. Signal Processing, 42(5):1202-1217, 1994], dans [G. MATHEW V. U. REDDY, Adaptive estimation of eigensubspace , IEEE Trans. Signal Processing, 43(2):401-411, 1995], dans [E. OJA, Neural networks, principal components and subspaces , Int. journal of neural systems, 1(1):61-68, 1989] et dans [G. Xu, H. ZHA, G. H. GOLUB T. KAILATH, Fast algorithms for updating signal subspaces , IEEE Trans. Circuits Syst., 41(8):537- 549, 1994]. En particulier, il est indiqué dans les publications [L. Xu, Least mean square error reconstruction principle for self-organizing neural nets , Neural Networks, 6:627-648, 1993] et [B. YANG, Projection Approximation Subspace Tracking , IEEE Trans. _ Signal Processing, 44(1):95-107, 1995], que la méthode classique de Oja divulguée dans [E. OJA, Neural networks, principal components and subspaces , Int. journal of neural systems, 1(1):61-68, 1989] peut être considérée comme une approximation d'un algorithme de descente de gradient d'une fonction d'erreur quadratique. Un certain nombre de méthodes de poursuite de sous-espace plus rapides ont été développées en combinant l'approche de descente de gradient avec l'hypothèse connue sous le nom de projection approximation , comme décrit dans [K. ABED-MERAIM, A. CHKEIF Y. HUA, Fast orthonormal PAST algorithm , IEEE Signal Proc. Letters, 7(3):60-62, 2000], dans [R. BADEAU, K. ABED-MERAIM, G. RICHARD B. DAVID, Sliding Window Orthonormal PAST Algorithm , Proc. of Int. Conf on Acoustic, Speech and Signal Processing, 5, 261-264. IEEE, 2003], dans [S. C. DOUGLAS, Numericallyrobust adaptive subspace tracking using Householder transformations , Proc. of Sensor Array and Multichannel Signal Proc. Workshop, 499-503. IEEE, 2000], dans [Y. MIAO Y. HUA, Fast subspace tracking and neural network learning by a novel information criterion , IEEE Trans. Signal Processing, 46(7):1967-1979, 1998], et dans [B. YANG, Projection Approximation Subspace Tracking , IEEE Trans. Signal Processing, 44(1) :95-107, 1995]. D'autres techniques connues reposent sur un moyennage des espaces signal et bruit, comme décrit dans [R. D. DEGROAT, Noniterative subspace tracking , IEEE Trans. Signal Processing, 40(3):571-577, 1992], sur le principe du maximum de vraisemblance comme décrit dans [C. RIOU T. CHONAVEL, Fast adaptive eigenvalue decomposition: a maximum likelihood approach , Proc. of Int. Conf on Acoust., Speech, and Signal Proc., 35653568, IEEE, 1997], sur l'analyse d'opérateurs restreints comme décrit dans [C. S. MACINNES, Fast, accurate subspace tracking using operator restriction analysis , Proc. of Int. Conf on Acoustic, Speech and Signal Processing, 1357- 1360. IEEE, 1998], ou sur la théorie des perturbations comme décrit dans [B. CHAMPAGNE Q.-G. Liu, Plane rotation-based EVD updating schemes for efficient subspace tracking , IEEE Trans. Signal Processing, 46(7):1886-1900, 1998]. Un résumé de techniques plus anciennes est proposé dans la publication [P. COMON G. H. GOLUB, Tracking a few extreme singular values and vectors in signal 10 processing , Proc. of IEEE, 78, 1327-1343, 1990]. Un autre algorithme de poursuite de l'espace signal, baptisé SP, a encore été proposé par C.E. Davila dans la publication [C. E. DAVILA, Efficient, high performance, subspace tracking for time-domain data , IEEE Trans. Signal Processing, 48(12):3307-3315, 2000]. Comme dans le cas de l'algorithme PAST divulgué dans [B. YANG, Projection Approximation Subspace Tracking , IEEE Trans. Signal Processing, 44(1):95-107, 1995], ou de l'algorithme NIC décrit dans [Y. MIA0 Y. HUA, Fast subspace tracking and neural network learning by a novel information criterion , IEEE Trans. Signal Processing, 46(7):19671979, 1998], l'estimation de l'espace signal est considérée comme un problème d'optimisation. Cependant, au lieu de déterminer de façon approximative le maximum global de la fonction à optimiser, l'algorithme SP détermine de façon exacte le maximum de cette même fonction sur un sous-espace de dimension limitée. Cet algorithme ne se prête qu'à l'analyse de séries temporelles, et présente un coût en calcul élevé. Sa complexité globale est O(nrz) (où n est la dimension des vecteurs de données observés et r est la dimension de l'espace signal), alors qu'un certain nombre d'algorithmes existants (c'est le cas par exemple de PAST, FST, NIC et OPAST) ne requièrent que O(nr) opérations à chaque instant. Toutes les publications citées sont incorporées dans la présente demande par référence. La présente invention vise notamment à proposer un procédé d'analyse spectrale à la fois performant et qui puisse être de coût en calcul moindre que l'algorithme SP, entre autres. Selon l'un de ses aspects s'appliquant notamment à l'analyse spectrale de données temporelles, l'invention a pour objet un procédé de poursuite d'un espace signal de dimension r inférieure à la dimension n des vecteurs de données, comportant l'étape consistant à calculer de manière récursive la matrice spectrale W(t) (en anglais subspace weighting matrix) de dimensions n x r engendrant ledit espace, par modification d'une matrice spectrale antérieure W(t z), où i peut être égal à 1. La matrice W(t) est avantageusement obtenue en calculant la solution exacte d'un problème d'optimisation contraint à un sous-espace de dimension comprise entre r+l et n. La matrice W(t) peut être mise à jour par addition d'une matrice de rang faible M(t), selon la récurrence W(t) = W(t r) + M(t). Dans un exemple de mise en oeuvre de l'invention, le procédé comporte le calcul de matrices intermédiaires W(t) et Z(t) et la mise à jour des matrices W(t) et Z(t). Le calcul de W(t) et Z(t) peut reposer sur un lemme d'inversion matricielle par ajout de lignes et de colonnes. La matrice W(t) peut être orthonormée à chaque itération, et la complexité algorithmique peut être de l'ordre de nr opérations. L'invention permet d'éviter de recourir à un algorithme de recherche d'une solution approchée dans un espace de dimension n comme c'est le cas dans le procédé PAST mentionné ci-dessus. De plus, dans un exemple de mise en oeuvre de l'invention, la mise à jour de W(t) peut s'effectuer par une récurrence simple, par ajout de la matrice M(t), sans multiplication matricielle complexe, comme c'est le cas dans l'algorithme SP. Le fait de rechercher la matrice M(t) dans un sous-espace de dimension limitée permet notamment de calculer plus rapidement la matrice W(t). La matrice de mise à jour M(t) peut être par exemple de rang 1, 2 ou 3, et de préférence de rang 1 si l'on cherche un suivi rapide de l'espace signal. Les données peuvent être des données audio, vidéo, ou d'imagerie, et plus généralement issues de mesures physiques (par exemple mécaniques). Dans le cas notamment où les données sont des données audio, celles-ci sont par exemple issues d'un signal audio éventuellement réhaussé, décomposé en sous-bandes, notamment par un banc de filtres en cosinus modulés, et dont le bruit a été préalablement blanchi par filtrage. Au moins l'une des sous-bandes peut être analysée par la mise en oeuvre du procédé d'analyse haute résolution selon l'invention. La décomposition en sous-bandes peut encore s'effectuer de toute autre manière appropriée. Le procédé selon l'invention peut trouver de nombreuses applications, notamment dans l'analyse spectrale, par exemple l'analyse automatique du rythme musical, le débruitage, notamment d'un signal de parole en vue d'un éventuel traitement ultérieur tel que la reconnaissance de la parole. Le procédé peut également s'appliquer au codage de données, notamment dans le cadre de la compression de données. La connaissance de la décomposition en sinusoïdes d'un signal donné peut permettre de le coder plus efficacement pour générer un fichier compressé. L'invention peut être appliquée à la séparation de la partie sinusoïdale et de la partie bruitée d'un signal, le bruit pouvant être éliminé par une projection sur l'espace signal. Cela peut être utile par exemple pour prétraiter la voix avant de mettre en oeuvre un procédé de reconnaissance vocale. L'invention peut également s'appliquer à des données non audio, par exemple des données vidéo ou d'imagerie médicale ou sismique, dans le but par exemple de restaurer une image. L'invention pourra être mieux comprise à la lecture de la description détaillée qui va suivre, d'exemples non limitatifs de mise en oeuvre de celle-ci, et à l'examen du dessin annexé, sur lequel: - la figure 1 est un schéma en blocs illustrant le pseudo-code d'un exemple de procédé d'analyse spectrale conforme à l'invention, - la figure 2 est un schéma en blocs représentant un exemple de pré-traitement du signal, et -la figure 3 est un schéma en blocs représentant un exemple de post- traitement du signal. 2888963 7 Analyse HR On suppose dans la suite que les données que l'on souhaite analyser se présentent sous la forme d'une séquence de vecteurs {x(t)}tEZ, de dimension n. Les données sont par exemple des données numériques résultant d'une conversion analogique/numérique puis d'un échantillonnage du signal à analyser. Le choix de la dimension n peut dépendre de la nature des données et de leur variabilité dans le temps. Les données sont par exemple des données audio, et peuvent avoir subi un pré-traitement, comme détaillé plus loin. Le procédé d'analyse haute résolution selon l'invention est adaptatif, c'est-à-dire que l'on tient compte des informations passées pour mettre à jour l'espace signal. Deux vecteurs successifs de la séquence {x(t)}tE z comportent une information redondante. La matrice formée par une séquence de vecteurs successifs peut ainsi présenter une forme prédéfinie, par exemple de Hankel ou de Toeplitz. Chaque vecteur est par exemple un segment de signal extrait d'une fenêtre temporelle de longueur n. La séquence de vecteurs x(t) peut correspondre à un ensemble de fenêtres décalées les unes des autres de un voire plusieurs échantillons. L'analyse de chaque fenêtre pourra ainsi reposer au moins partiellement sur l'analyse d'une fenêtre précédente. L'espace signal est défini comme l'espace propre principal de dimension r La matrice de corrélation est définie selon le type de fenêtre que l'on applique au signal (exponentielle ou tronquée). La fenêtre exponentielle permet d'atteindre une complexité plus faible, alors que la fenêtre tronquée offre un meilleur suivi des variations du signal. Dans ce qui suit, H désigne le conjugué hermitien. Fenêtre exponentielle: Cxx(t) _ E/3-ux(u)x(u)H, où 0 Cette matrice peut être mise à jour selon la récurrence: c. ) = /3 C (t 1)+x(t)x(t) H. (1) Fenêtre tronquée: Cxx(t) _ E /3t-x(u)x(u)H. (2) u=t-1+ 1 où 0 r est la longueur (finie) de la fenêtre. Le cas f = 1 correspond à une fenêtre rectangulaire, ou fenêtre glissante. Cette matrice peut être mise à 5 jour à l'aide de la récurrence: C (t) =C(t 1)+x(t)x(t)H J3'x(t l)x(t l)H. (3) Les équations (1) et (3) peuvent toutes deux être écrites sous la forme C (t) = J3 C(t 1)+x (t) J x(t)H. (4) où x (t) et J sont définis selon le type de fenêtre: -dans le cas de la fenêtre exponentielle: J= 1 (5) x(t) = x(t) (6) - dans le cas de la fenêtre tronquée: J= o QI (7) x(t) = [x(t) x(t -1)1 (8) Soit q le rang de la mise à jour impliquée dans l'équation (4). Puisque q =1 dans le cas de la fenêtre exponentielle et q = 2 dans le cas de la fenêtre tronquée, q caractérise la forme de la fenêtre. En particulier, x(t) est une matrice de dimensions n x q et J est une matrice de dimensions q x q. Il est connu qu'une matrice orthonormée W(t) de dimensions n x r engendre le sous-espace dominant de dimension r de la matrice C(t) si et seulement si elle maximise le critère J, défini comme V(W(t))= trace( W(t)''C (t) W(t))I Le maximum de ce critère est égal à la somme des r plus grandes valeurs 25 propres de C.,,(t). Cette maximisation n'est pas recherchée, conformément à un aspect de l'invention, parmi toutes les matrices orthonormées, car la complexité serait O(n2r) et ne conduirait pas à une récurrence simple entre W(t) et W(t-1). 15 La recherche s'effectue dans l'espace image de W(t-1), plus quelques directions de recherche additionnelles. L'espace image de la matrice W(t) (de dimension r) est obtenu comme un sous-espace de l'espace de dimension (r + p + q 1) (avec p = 1 ou 2) engendré par la matrice V(t) _ [W(t -1),x(t) (9) où la matrice x(t), de dimensions n x (p + q - 1), contient le vecteur x(t) dans tous les cas, le vecteur x(t - 1) dans le cas de la fenêtre tronquée où q = 2, et éventuellement le vecteur Cxx(t - 1) x(t) dans le cas p = 2. Soit W(t) une matrice orthonormée de dimensions n x (r + p + q - 1) 10 engendrant l'espace image de la matrice V(t) . Alors W(t) peut être écrite sous la forme W(t) = W(t)U(t) (10) où U(t) est une matrice orthonormée de dimensions (r + p + q - 1) x r. Dans ce cas J(W(t)) = trace(U(t)H -Ce(t)U(t)) où C, (t) est la matrice de dimensions (r + p + q - 1) x (r +p + q - 1) C,,, (t) = W (t)H Cy,(t)W (t). (12) Maximiser l'équation (11) revient à chercher la matrice U(t) qui, parmi toutes les matrices orthonormées U(t) de dimensions (r + p + q - 1) x r, engendre le sous-espace dominant de dimension r de la matrice C,, (t). Ainsi la matrice W(t) peut être suivie en calculant: - une base orthonormée W(t) de l'espace image de V(t) , - la matrice C,,, - (t) = W (OH Cxx (t)W - (t) , - une matrice orthonormée U(t) de dimensions (r + p + q - 1) x r engendrant le sous-espace dominant de dimension r de C,,,- (t) , - la matrice W(t)= W(t)U(t) . U(t) peut être obtenue comme une base du complémentaire orthogonal du sous- espace dominant de dimension (p + q 1) de la matrice Z(t) = (t)-' . Ainsi, dans un exemple de mise en uvre de l'invention, on calcule i(t) et son sous-espace dominant de dimension (p + q 1), et on calcule U(t) comme une base du complémentaire orthogonal de ce sous-espace. Comme cela sera montré plus loin, ce calcul peut être effectué en mettant à jour l'inverse Z(t) de la matrice de corrélation compressée Cyy(t), de dimensions r x r: Cy,(t)=W(t)HC (t)W(t). (13) Selon un exemple de mise en oeuvre de l'invention, W(t) peut être mise à jour selon un procédé comprenant les étapes de calcul de W- (t) , de Z(t) , la mise à jour de W(t) et la mise à jour de Z(t). Calcul de W(t) On considère la matrice y(t) = W(t -1)H-x(t) de dimensions r x (p + q 1), et on pose e(t) = x(t) W(t -1) y(t). (14) La matricee(t) , de dimensions n x (p + q 1), est orthogonale à l'espace image de W(t 1). Soit a- (t) une racine carrée de la matrice e(t)H e(t) , de dimension (p+q 1) x (p+q 1): 1 _ 6(t) _ (e(t)H e(t)) 2 = (x_(t)H x(t) _y(t)H _y(t)) z. (15) Ci-dessous, on suppose que a(t) est inversible, car même si o(t) se trouve être singulière, elle peut être rendue inversible en supprimant une ou plusieurs colonnes dans x(t) . La matrice W(t) = [W (t -1), e(t)6(t)-' j (16) est orthonormée. En particulier, V(t) peut être écrite sous la forme V(t) = W(t)R(t) (17) où R- (t) est une matrice inversible de dimensions (p + q 1) x (p + q 1) : 0(p+q-1) X r (18) Calcul de Z(t) Comme cela a été mentionné plus haut, la matrice Z(t) est définie comme l'inverse de la matrice ',,,(t) de dimensions (p + q 1) x (p + q 1) définie dans l'équation (12). Par conséquent, calculer Z(t) revient à calculer et inverser C,,,,(t). Calcul de C,,- (t) : en substituant l'équation (17) dans l'équation (12), on obtient C (t) = R(t)-H C' y,,(t)R(t)-' où C',,y(t) est une matrice de dimensions (r + p + q 1) x c' Y,(t) =V(t)HCxr(t)V(t)É Soient alors a(t) = x(t)H x(t) y(t) = w (t -1)H x(t) x'- (t) = C (t -1)x- (t) y'- (t) = w -1)H x'(t). En substituant les équations (4) et (9) dans l'équation (20), on obtient (21) (22) (23) (24) C'ri (t) = Cn, (t) ygt)H ygt) c (t) (25) (26) où Cy, (t) = fCn,(t 1)+ y(t)Jy(t)H y - (t) = /3 y (t) + y(t)J a- (t) (27) cy,(t) _ /3x(t)Hx'(t)+ a(t)H.1 a(t). (28) Inversion de C,,,, (t) : la matrice C,ry (t) peut être obtenue à partir de C yy(t) à l'aide de l'équation (19). Ainsi, inverser Cyy (t) nécessite d'inverser Cri (t). Cette dernière opération peut être effectuée en inversant d'abord le coin supérieur gauche de dimensions r x r de la matrice C,ry(t), noté Cyy (t). Inversion de C ,(t) : supposons que la matrice Cyy(t-1) soit inversible et posons Z(t-1) = Cyy(t-1)-'. Soit h(t) = Z(t -1)y(t). (29) Le lemme d'inversion matricielle [R. A. HORN, C. R. JOHNSON, Matrix analysis , Cambridge University Press, 1985, pp. 18-19] montre que Cn,(t) est inversible si et seulement si Pr' + y(t)H h(t) est inversible. Dans ce cas, soit y(t) = (/3 J-' +y(t)Hh(t))-1 Alors, i (t) C , (t) -' satisfait Z(t) = (Z(t -1) h(t)y(t)h(t)H). Inversion de C,(t) : inverser l'équation (19) implique Z(t) = R(t) Ce (t)-1 R(t)H. (32) Soit h(t) = Z(t) y"(t) . (33) Pour inverser la matrice C'y, (t) , on utilise le lemme suivant: Lemme d'inversion matricielle par ajout de lignes et de colonnes. Soient un entier r et une matrice inversible C de dimensions r x r. Soient un entier d, une matrice A de dimensions r x d, une matrice B de dimensions d x r, et une matrice D de dimensions d x d. On considère la matrice C = D, de dimensions (r + d) X (r + d). Alors C est inversible si et seulement si la matrice D BC-IA, de dimensions d x d, est inversible. Dans ce cas, on pose F = (D BC-14-1. Alors C- -C-' +C-'AFBC-' C-1AI' _ - rBC ' F (30) (31) En appliquant le lemme avec C = Cy,(t), A = y (t), B = y^(t)H, et D = (t) , on montre que C' (t) est inversible si et seulement si c, (t) y" (t)H h(t) est inversible. Dans ce cas, C'ri (t)-1 = Z(t) + h(t)y(t)h(t) H y(t)h(t) h(t)y(t)] y(t) H (34) où y(t) est la matrice définie positive de dimensions p x p y(t) _ (c,(t) yL(t)H h(t))-' É (35) En substituant les équations (18) et (34) dans l'équation (32), on obtient Z(t) = Z'(t) - g(t)H g(t) y'(t) (36) où h'(t) = h(t) y(t) (37) Z' (t) = Z(t) + h' (t)y(t)h' (t)H (38) g(t) = h' (t)y(t)6(t)H (39) y'(t) = c (t)y(t)6(t)H. (40) Mise à jour de W(t) Soit (t) une matrice orthonormée de dimensions (r + p + q 1) x (p + q 1), dont les colonnes engendrent le sous-espace dominant de dimension (p + q 1) de la matrice Z(t), de dimensions (r +p + q 1) x (r + p + q 1). Il existe une matrice À(t) définie positive de dimensions (p + q 1) x 20 (p + q 1), dont les valeurs propres sont égales aux p + q 1 plus grandes valeurs propres de Z(t) , telle que Z(t)O(t) = O(t)Z(t) (41) Soit t) la matrice de dimensions r x (p + q 1) contenant les r premières colonnes de 0(t), et z(t) la matrice de dimensions (p + q -1) x (p + q 1) contenant ses 25 p + q 1 dernières lignes: 0(t)T = [9(t)T, z(t)T]É (42) La factorisation polaire de z(t) est définie comme z(t) = p(t)e(t) (43) où p(t) est définie positive et 0(t) est orthonormée. Soit f (t) = tp(t)6(t)H (44) f'(t) = f(t)(I + p(t))-' . (45) Alors, on peut vérifier que la matrice de dimensions (r + p + q 1) x r U(t) = Ir - f (t)f(t)H _f(t)H (46) est orthonormée et satisfait U(t) H O(t) = 0. Par conséquent, U(t) est une base orthonormée du sous-espace mineur de dimension r de la matrice Z(t) . En substituant les équations (16) et (46) dans l'équation (10), on obtient une récurrence pour la matrice W(t) : W(t) = W(t -1) e (t)f(t)H (47) où e (t) = e(t)6(t)-' + W(t -1) f (t) . Finalement, en substituant l'équation (14) dans cette dernière définition, on 15 obtient e (t) = x(t)6(t)-' W(t -1) y (t) (48) Yo(t) =Y(t)6(t)-' -7 (t) É (49) Mise à jour de Z(t) En substituant les équations (10) et (12) dans l'équation (13), on obtient Cam, (t) = U(t) H Cam, (t)U(t) . Comme la matrice orthonormée U(t) engendre un sous-espace invariant de Z (t) , elle engendre aussi un sous-espace invariant de son inverse Cam, (t) . Par conséquent, 25 cette dernière équation implique Z(t) = U(t)H Z(t)U(t) . (50) où 2888963 15 En substituant les équations (45), (44) et (42) dans l'équation (46), on montre U(t) f (t) H - (D(t)b (t) H f ' (t) H. (51) Ainsi, en substituant les équations (51), (36), (41) dans l'équation (50), on obtient Z(t) = Z (t) + g (t).f (t)H + f (t)8(t)H (52) où g t) = g(t) + .f (t)(Y (t) 0 (t)X (t)e (t)H) . (53) Implémentation Un exemple de pseudo-code du procédé d'analyse spectrale selon l'invention est représenté à la figure 1. On construit à l'étape 10 le vecteur de données x(t), t étant l'indice d'itération, puis le vecteur de données compressées y(t) à l'étape 11, la direction de recherche x'(t) à l'étape 12 et le vecteur de données compressées y'(t) de x'(t)à l'étape 13. On détermine à l'étape 14 si q = 2 et dans l'affirmative on effectue les étapes 15 à 19, spécifiques à la fenêtre tronquée. Dans la négative, on va à l'étape 20. A l'étape 22, on détermine si p = 2 et dans l'affirmative on exécute les étapes 23 à 27, destinées à affiner l'estimation de l'espace signal. Dans la négative, on va à l'étape 28 de calcul de l'écart type de l'erreur de prédiction 6(t), suivi des étapes 29 à 51. L'étape 42 fait intervenir le calcul des p+q-1 principaux vecteurs et valeurs propres de Z(t). L'étape 44 correspond à la décomposition polaire de z(t). L'étape 48 correspond notamment au calcul du vecteur de mise à jour e'(t) de l'espace signal. L'étape 52 force la symétrie de Z(t), de manière à tenir compte des erreurs d'arrondi pouvant survenir lors du calcul. que 2888963 16 Après l'étape 52, l'indice d'itération t est augmenté d'une unité et l'on retourne à l'étape 10. La complexité de l'algorithme ainsi obtenu peut être évaluée en calculant le nombre de Multiplications / Accumulations (MACs) qu'il comporte. Dans le cas de la fenêtre exponentielle, son coût est ainsi de 4nr MACs si p = 1, ce qui est inférieur ou égal à celui de FST, NIC et OPAST, et 7nr si p = 2, ce qui reste un ordre de grandeur plus faible que celui des algorithmes de Karasalo et LORAF ou SP. Dans le cas de la fenêtre tronquée, son coût global est 8nr MACs si p = 1, ce qui est généralement plus faible que celui de TW-FAPI, et 11 nr si p = 2, ce qui reste plus 10 faible que celui de SHSVD2, SWASVD2 et SW-OPAST. Le coût du calcul des vecteurs x'(t) = Cxx(t 1) x(t), x"(t) = Cxx(t 1) 2 x(t) et z t I) = Cxx (t 1) x(t 1) peut être réduit de 0(n2) à 0(n) en utilisant la technique décrite dans [C.E. DAVILA Efficient, highperformance, subspace traking for time-domain data , IEEE Trans. Signal Processing, 48(12) :3307-3315, décembre 2000], qui exploite la structure Hankel de la matrice de données. Le procédé peut être modifié pour suivre le sous-espace mineur de la matrice de corrélation, auquel cas la matrice U(t) doit engendrer le sousespace mineur de dimension r de C (t) , au lieu de i(t) . Système d'analyse de données audio On a représenté à la figure 2 un exemple de système d'analyse d'un signal d'entrée audio. Le signal entrant est partitionné en segments adjacents (déterminés de façon à ce que l'ordre de modélisation puisse être choisi constant sur la durée de chaque segment), puis filtré de façon à rehausser les hautes fréquences par rapport aux basses fréquences, et découpé ensuite en sous-bandes à l'aide d'un banc de filtres. Le réhaussement des hautes fréquences peut être effectué par le biais d'un filtre de pré-accentuation, par exemple de fonction de transfert H(z) =1-0.98z"'. Lors de la décomposition en sous-bandes, afin d'obtenir une bonne estimation des paramètres fréquentiels dans chaque sous-bande, il est souhaitable d'utiliser des filtres d'analyse suffisamment réjecteurs pour que la puissance du signal dans la bande atténuée ne dépasse jamais le niveau de bruit dans la bande passante. De préférence, la largeur de bande en hautes fréquences n'excède pas 2500 Hz, et la longueur des filtres d'analyse ne dépasse pas la durée moyenne de stationnarité du signal, par exemple de l'ordre de 25 ms. La technique de découpage peut faire intervenir un banc de filtres en cosinus modulés (CMF). Ce type de banc de filtres a déjà été utilisé en pré-traitement de méthodes sous-espace dans [K. HERMUS, W. VERHELST ET P. WAMBACQ Psycho acoustic modeling of audio with exponentially damped sinusoids , Proc. of ICASSP '02, Volume 2, pages 1821-1824, IEEE, 2002], aboutissant à un découpage uniforme du spectre. Cependant, il est possible de construire un banc de filtres non uniforme en regroupant des voies adjacentes (c'est-à-dire en les sommant), comme cela a été proposé dans [J-J. LEE ET B.G. LEE, A design of non uniform cosine modulated filter banks , IEEE Trans. Circuits Syst. II, 42(11):732-737, novembre 1995]. Chaque sous-bande ainsi obtenue peut alors être décimée d'un facteur égal au nombre total de voies du banc CMF, divisé par le nombre de voies fusionnées. Par exemple, dans le cas d'un banc CMF contenant 32 voies, une sous- bande obtenue en fusionnant 4 voies adjacentes peut être décimée d'un facteur 8. Sous certaines conditions, le banc de filtres ainsi obtenu satisfait trois propriétés: les filtres d'analyse et de synthèse sont passe-tout dans leur bande passante ^ le banc de filtres d'analyse/synthèse vérifie la condition d'annulation de repliement ^ la fonction de distorsion est identique à celle du banc CMF. De plus, après décimation critique, tous les filtres d'analyse deviennent passe-tout. Ainsi, si le bruit dans le signal en entrée du banc de filtres est blanc, le bruit dans les signaux de sous-bande est également blanc, ce qui constitue une propriété souhaitable dans le cadre de l'analyse HR. La table ci-dessous présente un exemple de banc de filtres obtenu par cette méthode. Le signal en entrée est échantillonné à 44100 Hz et l'on se base sur un banc CMF contenant M=32 voies, numérotées de 0 à M-1. Après regroupement des voies, le banc de filtres contient alors 8 sous- bandes non-uniformes. En pratique, la Sème sous-bande n'est pas traitée par l'analyse HR, car l'on suppose dans cet exemple qu'il n'y a pas ou peu de sinusoïdes au-delà de 11025 Hz. Sous-bande 0 1 2 3 4 5 6 7 Voies 0 1 2-3 4-5 6-7 8-11 12-15 16-31 fusionnées Intervalle 0-690 690-1380 1380-2755 2755- 4135- 5515- 827011025- (Hz) 4135 5515 8270 11025 22050 Décimation 32 32 16 16 16 8 8 2 L'étape suivante dans la figure 2 consiste à blanchir le bruit additif contenu dans les signaux de sous-bandes. Pour cela, la densité spectrale de puissance du bruit peut être estimée pour en déduire un filtre blanchisseur dans chaque sous-bande. Une méthode d'estimation de la densité du bruit en présence de sinusoïdes a été proposée dans [R. VON SACHS, Peak-insensitive non-parametric spectrum estimation , Journal of time series analysis, 15: 429-452, 1994] , ce qui comprend le calcul d'un périodogramme lissé à l'aide d'un noyau insensible aux pics fréquentiels. Ce périodogramme est obtenu comme la solution d'une équation implicite qui peut être résolue de manière itérative. Sous certaines hypothèses, il est démontré que cet estimateur est asymptotiquement consistant. L'implémentation de cette technique itérative s'avère cependant assez coûteuse en temps de calcul. De préférence, le périodogramme est lissé non pas par filtrage linéaire, mais à l'aide d'un filtre de rang, que l'on peut choisir insensible aux pics fréquentiels. Ainsi, le résultat du filtrage sera comparable, mais la résolution itérative d'une équation implicite aura été évitée. Une fois le périodogramme estimé dans chaque sous-bande, il est possible d'en déduire la fonction d'autocovariance du bruit. Des méthodes classiques de prédiction linéaire [L. L. SCHARF, Statistical signal processing , Addison Wesley, New York, 1991] permettent alors de déterminer dans chaque sous-bande un filtre à Réponse Impulsionnelle Finie (RIF) dont l'action est de blanchir le bruit. Les filtres ainsi obtenus sont finalement appliqués aux signaux de sous-bandes. Traitements applicables dans les sous-bandes Des traitements peuvent être effectués directement sur les signaux de sous- bandes. Par exemple, dans le cas des applications de codage, il suffit de coder les paramètres des exponentielles complexes (ou des sinusoïdes) estimées dans les sousbandes (la partie bruitée pouvant par exemple être modélisée par un processus tel que décrit par [M. KAHRS et K. BRANDENBOURG dans Applications of Digital Signal Processing to Audio and Acoustics , Kluwer Academic Press, Dortrecht, Pays-Bas, 1998]). Les signaux de sousbandes peuvent ensuite être resynthétisés à partir de ces paramètres, puis injectés dans le système de reconstruction du signal de la figure 3. Il est également possible d'envisager des opérations de débruitage et d'extraction du bruit, qui sont complémentaires l'une de l'autre: il s'agit de séparer la partie sinusoïdale du signal de sa partie bruitée. L'extraction de la partie bruitée du signal peut avoir plusieurs applications, dont l'estimation du rythme musical. Débruitage du signal et extraction du bruit La paramétrisation des signaux de sous-bandes permet de séparer leur partie sinusoïdale de leur partie bruitée. En effet, dans chaque fenêtre d'analyse, il est possible de resynthétiser la partie sinusoïdale à partir des paramètres estimés Lia partie bruitée est alors constituée du résiduel, obtenu en soustrayant ce signal de synthèse au signal original. Une simple technique d'addition-recouvrement permet ensuite de reconstruire les signaux séparés à partir des segments obtenus sur chaque fenêtre. Il est possible d'utiliser par exemple des fenêtres de Hanning avec un recouvrement égal à la moitié ou au quart de la longueur des fenêtres. Les signaux séparés peuvent être injectés dans le système de reconstruction de la figure 3. Le procédé d'analyse HR et les pré ou post traitements éventuels associés peuvent être mis en oeuvre dans un ordinateur programmé pour effectuer le traitement des données conformément à l'invention. Ce traitement peut s'effectuer en temps réel ou en temps différé. Le dispositif servant à mettre en oeuvre l'invention comporte avantageusement une mémoire permettant de stocker les échantillons de signal à analyser et les données intermédiaires résultant de l'exécution du procédé d'analyse HR. Le procédé d'analyse peut encore être mise en oeuvre par un circuit 30 électronique spécialisé, comportant par exemple au moins un micro contraeur ou micro processeur et les composants périphériques correspondants. L'invention peut encore être mise en oeuvre à distance, par le biais d'un réseau informatique et/ou de téléphonie. Dans ce cas, les données à traiter sont par exemple envoyées à un serveur distant par le réseau, ce serveur effectuant l'analyse HR, puis les paramètres du signal peuvent être retournés, par le même réseau le cas échéant. L'invention peut s'appliquer à l'analyse du rythme musical, à la reconnaissance d'instruments et/ou de notes de musique, par exemple en vue de la transcription automatique de la musique. L'analyse HR de données audio peut être utile par exemple pour l'indexation automatique du contenu musical diffusé par une station, par exemple de radio ou de télévision, en vue par exemple du calcul des droits à verser aux auteurs dont les oeuvres sont diffusées. L'analyse HR de la voix peut être utile afin par exemple de débruiter le signal pour obtenir une meilleur écoute, ou d'améliorer un traitement ultérieur tel que la reconnaissance vocale, ou encore de déterminer des temps de parole, ce qui peut servir notamment à contraer le bon déroulement d'une campagne électorale. Bien que les exemples ci-dessus concernent avantageusement des données audio, l'invention s'applique encore à des données non audio, par exemple vidéo. L'invention peut par exemple être appliquée au traitement de données issues de l'imagerie médicale, par exemple de la RMN, ou plus généralement issues de mesures physiques (par exemple mécaniques ou sismiques). REVENDICATIONS 1. Procédé de poursuite d'un espace signal de dimension r inférieure à la dimension n des vecteurs de données, comportant l'étape consistant à calculer de manière récursive la matrice spectrale W(t) de dimensions n x r engendrant ledit espace, par modification d'une matrice spectrale antérieure W(t z), où r peut être égal à 1. 2. Procédé selon la revendication 1, caractérisé par le fait que la matrice W(t) est obtenue en calculant la solution exacte d'un problème d'optimisation contraint à un sous-espace de dimension comprise entre r+l et n. 3. Procédé selon l'une des revendications précédentes, caractérisé par le fait que la matrice W(t) est mise à jour par addition d'une matrice de rang faible M(t), notamment égal à 1,2 ou 3, selon la récurrence W(t) = W(t r) + M(t). 4. Procédé selon l'une quelconque des revendications précédentes, caractérisé par l'utilisation des matrices intermédiaires W(t) et Z(t) dont le calcul repose sur un lemme d'inversion matricielle par ajout de lignes et de colonnes et par la mise à jour des matrices W(t) et Z(t). 5. Procédé selon l'une quelconque des revendications précédentes, caractérisé par le fait que la matrice W(t) est orthonormée à chaque itération. 6. Procédé selon l'une quelconque des revendications précédentes, caractérisé par le fait que les données sont des données audio, vidéo, ou d'imagerie, et plus généralement issues de mesures physiques, notamment mécaniques. 7. Procédé selon l'une quelconque des revendications précédentes, caractérisé par le fait que les données sont issues d'un signal audio éventuellement réhaussé, décomposé en sous-bandes, notamment par un banc de filtres en cosinus modulés, et dont le bruit a été préalablement blanchi par filtrage. 8. Procédé selon l'une quelconque des revendications précédentes, caractérisé par le fait qu'il est appliqué à la séparation de la partie sinusoïdale et de la partie bruitée d'un signal, ou au débruitage, notamment d'un signal de parole en vue d'en éventuel traitement ultérieur tel que la reconnaissance de la parole. 9. Procédé selon l'une quelconque des revendications précédentes, caractérisé par le fait qu'il est appliqué à l'analyse automatique d'un morceau de musique, notamment en vue de déterminer le rythme et/ou de le transcrire. 10. Procédé selon l'une quelconque des revendications précédentes, caractérisé par le fait qu'il est appliqué au codage de données, notamment d'un signal audio. 11. Procédé selon l'une quelconque des revendications précédentes, caractérisé par le fait qu'il est appliqué à l'analyse spectrale d'un signal, notamment un signal audio.