FR 2905489 A1 20080307 FR 0653507 A 20060829 1 Procédé d'estimation de phase pour la modélisation sinusoïdale d'un signal numérique 5 La présente invention se situe de manière générale dans le domaine du traitement du signal, et plus précisément dans le domaine de l'analyse, de la compression et de l'indexation des signaux. 10 L'invention concerne l'estimation de phase d'une composante sinusoïdale d'un signal numérique, basée sur l'analyse de la transformée de Fourier à court terme de ce signal. En modélisation sinusoïdale des signaux audionumériques, on représente un signal sonore par un ensemble d'oscillateurs dont les 15 paramètres d'amplitude et de fréquence varient lentement au cours du temps. Grâce à une analyse à court terme du signal dans les domaines temporels et fréquentiels, on détermine les valeurs de ces paramètres. Cette analyse est communément effectuée par une transformée de Fourier discrète à court terme du signal, c'est-à-dire calculée sur un nombre fini d'échantillons du 20 signal. Les pics du spectre en fréquence obtenu sont isolés, afin d'attribuer à chaque composante sinusoïdale correspondante, les valeurs d'amplitude, de fréquence et de phase associées. En ce qui concerne l'estimation de la phase d'un pic du spectre en fréquence du signal, la méthode la plus simple et la plus utilisée consiste à 25 calculer son argument. Cette méthode fonctionne correctement pour les composantes sinusoïdales stationnaires en fréquence, c'est-à-dire pour les composantes dont la fréquence est constante pendant toute la durée d'une fenêtre d'analyse. En revanche, pour une composante de fréquence variable pendant 30 cette durée, la phase estimée présente un biais. Pour des variations de fréquence marquées par rapport à la fenêtre d'analyse, le biais obtenu dans 2905489 2 l'estimation de la phase d'un pic est potentiellement important, entraînant une modélisation approximative du signal à analyser. Cette modélisation approximative est nuisible, par exemple dans un procédé de re-synthèse du signal, le signal re-synthétisé à partir de cette modélisation présente une mauvaise qualité auditive, s'il s'agit d'un signal sonore. Une précédente demande de brevet FR0553891 déposée par les demandeurs propose une solution pour compenser le biais obtenu dans une première estimation de phase d'un pic par la méthode précitée de calcul de l'argument du pic. Cette solution consiste à corriger la première estimation de phase, par un polynôme de degré au moins un en variation de fréquence estimée sur la fenêtre d'analyse du signal à analyser. Plusieurs modes de réalisation de cette solution sont possibles: Le premier mode utilise uniquement une Transformée de Fourier Discrète (TFD) du signal. Cependant les fréquences du spectre issu de la TFD ne correspondent pas aux vraies fréquences des composantes sinusoïdales du signal. Ceci a pour conséquence que l'estimation de phase corrigée par ce premier mode de réalisation présente parfois encore un biais sensible. Cette estimation des phases des composantes sinusoïdales du signal est donc sous-optimale, ce qui entraîne une modélisation perfectible du signal. Les autres modes utilisent une Transformée de Fourier à Temps Discret (TFTD), ce qui conduit à une estimation de phase corrigée permettant d'avoir une bonne modélisation du signal à analyser. Néanmoins dans ces autres modes, il est nécessaire d'appliquer une TFTD par composante sinusoïdale du signal, ce qui nécessite un temps de calcul plus important que dans le premier mode de réalisation. Or le temps de calcul est un facteur essentiel qui peut être préjudiciable dans certaines applications en temps réel, par exemple pour du codage audionumérique de type "sinusoïdal plus bruit". 2905489 La présente invention a pour but d'améliorer la technique antérieure en fournissant un procédé et un dispositif d'estimation améliorée de la phase d'une composante sinusoïdale d'un signal numérique non stationnaire. Dans 5 ce procédé et ce dispositif, on corrige la première estimation de phase obtenue par la méthode de calcul de l'argument du pic, par une valeur de biais tenant compte de l'écart de fréquence entre la fréquence estimée de la composante sinusoïdale et la fréquence discrète du spectre de Fourier correspondante. Le procédé selon l'invention a été testé en calculant cette 10 valeur de biais à partir d'une grille de valeurs prédéterminées, cette grille ayant été construite par avance pour une fenêtre d'analyse donnée et une fréquence d'échantillonnage donnée. Les résultats obtenus montrent une nette amélioration de l'estimation de phase par rapport aux techniques existantes utilisant uniquement la Transformée de Fourier Discrète. 15 A cette fin, l'invention propose un procédé d'estimation de phase d'une composante sinusoïdale d'un signal numérique non stationnaire, comportant les étapes de: - segmentation dudit signal en trames temporelles, 20 Transformée de Fourier Discrète d'une desdites trames, pondérée par une fenêtre d'analyse afin d'obtenir une valeur d'un pic du spectre en fréquence dudit signal correspondant à ladite composante sinusoïdale, estimation de fréquence sur ladite trame de ladite composante sinusoïdale, 25 première estimation de phase de ladite composante sinusoïdale en utilisant l'argument de ladite valeur de pic ainsi obtenue, caractérisé en ce qu'il comprend les étapes supplémentaires de : - estimation d'une valeur de biais en fonction de la différence de fréquence entre la fréquence précédemment obtenue à l'étape 30 d'estimation et la fréquence discrète dudit pic, et en fonction d'une 2905489 4 variation de fréquence estimée de ladite composante sinusoïdale entre le début et la fin de ladite fenêtre d'analyse, correction de ladite première estimation de phase par ladite valeur de biais. 5 Ce procédé permet d'obtenir une estimation de phase plus précise par rapport aux techniques connues utilisant uniquement la TFD comme moyen d'analyse spectrale, et plus rapide par rapport aux autres techniques, utilisant par exemple la TFTD. L'amélioration de la précision est d'autant plus notable que le signal est bruité. Par rapport aux applications nécessitant une analyse 10 de signal rapide et n'utilisant que la TFD, les applications utilisant le procédé selon l'invention se trouvent améliorées. Ainsi : Une application de suivi de partiel utilisant l'invention est améliorée car le contrôle du verrouillage de la phase entre deux analyses consécutives est meilleur, 15 Une application de re-synthèse d'un signal audionumérique à partir de ses partiels utilisant l'invention présente une qualité auditive améliorée car le signal re-synthétisé est plus proche de l'original, Dans une application de codage de signaux en signaux transitoires, sinusoïdaux et bruités utilisant l'invention, une meilleure estimation de 20 la phase des partiels conduit à une meilleure estimation du signal résiduel, ou de bruit, par soustraction du signal original et des composantes sinusoïdales estimées. Cela conduit à un signal résiduel de variance plus faible donc perceptivement plus facile à coder. Selon une caractéristique préférée, l'étape d'estimation d'une valeur de 25 biais utilise une grille de valeurs prédéterminées comprenant, pour des valeurs de différences de fréquences et de variations de fréquence prédéterminées, des valeurs de biais 2905489 B(8[u], A[v]) = ù arct oû B(8[u],[v]) est une valeur de biais de la grille de valeurs prédéterminées, 8[u] est une valeur de différence de fréquences prédéterminée, 5 A[v] est une valeur de variation de fréquence prédéterminée, n est l'indice d'échantillonnage de la trame, N est la longueur de la trame, Te est la période d'échantillonnage du signal, et w(n) est la valeur de la fenêtre d'analyse à l'échantillon n. 10 La grille utilisée pour calculer le biais par rapport à la première estimation de phase contient ainsi des valeurs très proches des valeurs théoriques de biais, pour une différence de fréquence prédéterminée entre la fréquence d'une composante sinusoïdale d'un signal et la fréquence correspondante dans le spectre de la TFD de ce signal, et pour une variation 15 de fréquence prédéterminée dans la fenêtre d'analyse du signal. Une telle grille prend un temps de calcul important à remplir mais n'est remplie qu'une seule fois, pour une fenêtre d'analyse donnée et une fréquence d'échantillonnage donnée. Selon une caractéristique préférée, dans l'étape d'estimation d'une 20 valeur de biais, ladite valeur de biais est obtenue par interpolation de valeurs prédéterminées de ladite grille. La valeur de biais obtenue à l'issue de l'étape d'estimation d'une valeur de biais est ainsi très proche de sa valeur théorique, ce qui donne une estimation de phase corrigée très précise. 2905489 6 L'invention concerne aussi un procédé de modélisation harmonique d'un signal numérique non stationnaire, utiEisant le procédé d'estimation de phase selon l'invention. L'invention concerne également un dispositif d'estimation de phase 5 d'une composante sinusoïdale d'un signal numérique non stationnaire, caractérisé en ce qu'il comporte des moyens de mise en oeuvre du procédé selon l'invention. L'invention concerne encore un programme d'ordinateur caractérisé en ce qu'il comporte des instructions de mise en oeuvre du procédé selon 10 l'invention, lorsque ledit programme est exécuté dans un système informatique. Les avantages de ce dispositif ou de ce programme d'ordinateur sont identiques à ceux mentionnés plus haut en relation avec le procédé selon l'invention. 15 D'autres caractéristiques et avantages apparaîtront à la lecture d'un mode de réalisation préféré décrit en référence aux figures dans lesquelles: - la figure 1 est un organigramme représentant les étapes du procédé selon un mode de réalisation de l'invention, 20 - la figure 2 est la représentation d'un signal numérique dont on cherche à modéliser les composantes sinusoïdales selon l'invention, - la figure 3 est la représentation de ce signal numérique sur une trame, - la figure 4 est un tracé permettant de comparer le procédé d'estimation de phase selon l'invention avec d'autres procédé d'estimation de phase de l'art 25 antérieur utilisant uniquement la Transformée de Fourier Discrète comme moyen d'analyse spectrale. Selon un mode de réalisation de l'invention, le procédé d'estimation de phase selon l'invention est utilisé dans un procédé plus général de 30 modélisation harmonique d'un signal numérique utilisant une analyse en 2905489 7 temps et en fréquence à court terme. Le schéma de principe de ce procédé est représenté à la figure 1, Les étapes al, a2, a3 et a5 correspondent au procédé de modélisation harmonique tel qu'il est habituellement implémenté, l'amélioration de 5 l'estimation de la phase du signal numérique selon l'invention étant implémentée à l'étape a4. Typiquement le signal numérique analysé est un signal sonore, que l'on modélise par un ensemble d'oscillateurs dont les paramètres d'amplitude et de fréquence varient lentement au cours du temps. Cependant le procédé selon 10 l'invention est applicable à tout signal numérique présentant les mêmes caractéristiques que ce type de signal sonore. Notamment les fréquences des composantes sinusoïdales du signal doivent être suffisamment espacées les unes des autres, afin de permettre son analyse spectrale. De plus les étapes al à a5 sont typiquement implémentées chacune 15 dans un module fonctionnel, à l'intérieur d'un équipement d'analyse de signaux numériques, par exemple un ordinateur. Ces modules fonctionnels utilisent des moyens logiciels ou des moyens matériels électroniques. Dans la première étape al, le signal numérique x(n) à modéliser, de 20 fréquence d'échantillonnage Fe, et d'indice d'échantillon n, est découpé en trames de quelques dizaines de millisecondes. Par exemple on travaille sur des trames de 512 échantillons pour un signal échantillonné à 16kHz (kilo Hertz). Comme représenté à la figure 2, l'analyse des composantes sinusoïdales du signal numérique se fait sur chacune des trames consécutives 25 du signal x(n), qui se recouvrent en partie car on choisit en général le pas d'avancement M des trames d'une valeur inférieure à la longueur N des trames. Ceci permet d'avoir des estimations des paramètres sinusoïdaux du signal de manière rapprochée au cours du temps. Il est à noter que par un souci de simplification de la figure 2, les trames d'indices p et p+l qui y sont 30 représentées présentent très peu d'échantillons. 2905489 8 On obtient ainsi un signal xp(n), représenté à la figure 3, où p est l`indice de la trame considérée, défini par la formule: xp(n) = x(n + p~) pour n E [0,N -1] où 5 p est l'indice de la trame considérée, M est le pas d'avancement des trames, N est la longueur de chaque trame, x(n) est la valeur du signal numérique à l'échantillon n, xp(n) est la valeur du signal numérique à l'échantillon n de la trame d'indice p. 10 Les étapes al à a5 sont dans cet exemple de réalisation effectuées de manière continue sur le signal x(n), c'est à dire sur les trames successives du signal x(n) pour les étapes a2 à a5. Dans la deuxième étape a2, on effectue une analyse spectrale du signal à court terme pour chaque trame du signal x(n). Typiquement on 15 applique au signal numérique fenêtré par une fenêtre d'analyse, une Transformée de Fourier Discrète. On obtient alors pour chaque trame d'indice p le spectre Xp(k), k étant un indice fréquentiel variant de 0 à N-1, déterminé par l'équation suivante: ,v 2 N X (k) _ x (n + ù)w(n + -) exp(ù j2~tnk) , p N v 2 2 N 2 20 où w(n) est une fenêtre d'analyse qui détermine la résolution fréquentielle de la transformée de Fourier. Cette fenêtre d'analyse est par exemple une fenêtre de type Blackman ou Hann, la seule condition requise sur cette fenêtre de longueur N, c'est-à-dire de même longueur qu'une trame de signal x(n), étant qu'elle soit symétrique en N/2, selon les équations suivantes: 25 w(N/2-n)=w(N/2+n), pour n variant de 1 à (N/2 - w(0)=0 Le choix particulier de cette fenêtre d'analyse et la permutation circulaire d'ordre N/2 effectuée sur le produit de xp(n) avec w(n), aussi appelée opération de "zero-phasing", avant la Transformée de Fourier Discrète, permet 2905489 9 à l'étape a3 d'estimer la phase des pics ou maxima locaux du spectre Xp(k) au milieu de l'intervalle d'analyse, c'est-à-dire à l'échantillon N/2 de la trame d'indice p dans le domaine temporel. L'intervalle d'analyse correspond en effet à la fenêtre d'analyse, il en a donc la même longueur, qui est aussi la taille de 5 la TFD. Dans la troisième étape a3, les pics, ou maxima locaux du spectre Xp(k), sont extraits. Seuls les pics les plus pertinents sont retenus, un critère étant par exemple de ne garder que les pics les plus énergétiques. A chacun des pics ainsi isolés sont ensuite affectées les valeurs en amplitude, en 10 fréquence et en phase correspondantes. En général, les fréquences des composantes sinusoïdales du signal à analyser ne coïncident pas avec les fréquences de la TFD utilisée qui, par définition, sont discrètes et liées à la fréquence d'échantillonnage Fe et à la taille N de la TFD par la relation : 15 Fk = F N où k est l'indice de fréquence du spectre Xp(k) et Fk la fréquence associée en Hertz. C'est pourquoi, la fréquence de chaque pic est re-estimée après-coup, ce qui permet aussi de re-évaluer l'amplitude du pic en fonction, par exemple, 20 de la fréquence du maximum local associé au pic et de la fréquence reestimée. De nombreuses méthodes existent pour re-estimer la fréquence du pic, par exemple les méthodes d'interpolation ou de réassignement à partir des raies TFD au voisinage du maximum local. Notamment la méthode du vocodeur de phase, divulguée dans la publication "Phase-locked vocoder" de 25 Miller Puckette en 1995 dans le cadre d'un atelier "Acoustics, Speech and Signal Processing" (ASSP) de l'organisme "Institute of Electrical and Electronic Engineer" (IEEE), donne une bonne estimation de la fréquence sur l'intervalle ou fenêtre d'analyse. Il en est de même avec la méthode de la dérivée, explicitée dans le document de thèse de décembre 2000 de Sylvain 30 Marchand, intitulée "Modélisation informatique du son musical (analyse, 2905489 10 transformation, synthèse)", de l'Université Bordeaux 1. Ainsi pour un pic d'indice fréquentiel K dans le spectre Xp(k on obtient une bonne estimation fK(N/2)de la fréquence réelle de la composante sinusoïdale correspondante du signal xp(n), au milieu de l'intervalle d'analyse, c'est-à-dire à l'échantillon N12 de la trame d'indice p. Dans la suite, on ne s'intéresse pour simplifier qu'à la composante sinusoïdale correspondant au pic d'indice fréquentiel K dans le spectre Xp(k), dont la fréquence estimée est fK (N/2) . Une estimation de l'amplitude du pic d'indice fréquentiel K est par exemple obtenue par la formule: 10 AK = où Xp(K W(FK - fK (Nt 2)) AK est l'estimation de l'amplitude du pic d'indice fréquentiel K, Xp(K) est la valeur de la raie du spectre Xp(k) à l'indice fréquentiel K, désigne la norme de la variable complexe z, 15 W(f) est la valeur de la transformée de Fourier de la fenêtre d'analyse w(n) à la variable f de fréquence, FK est la fréquence correspondant à l'indice fréquentiel K dans le spectre Xp(k), et fK(N/2)est la fréquence estimée à l'échantillon N/2 du pic d'indice 20 fréquentiel K dans le spectre Xp(k). En ce qui concerne l'estimation de la phase de la composante sinusoïdale de fréquence estimée siK (N/2), l'opération de "zero-phasing" effectuée à l'étape a2 permet de s'affranchir des rotations de phase dues aux écarts entre la fréquence réelle du pic à extraire et la fréquence de la raie TFD 25 qui lui est associée. Le document de thèse de Sylvain Marchand mentionné plus haut montre que la phase obtenue en prenant l'argument du pic est alors phase de la composante sinusoïdale correspondante au milieu de 2905489 11 l'intervalle d'analyse. On obtient ainsi pour le pic d'indice fréquentiel K, la phase estimée de la composante sinusoïdale correspondante à l'échantillon N/2 de la trame d'indice p, par la formule: Op (N/2) = Arg Xp(K)) 5 où K çop (Ni 2) est la phase estimée de la composante sinusoïdale considérée à l'échantillon N/2, et K est l'indice fréquentiel du maximum local correspondant dans le spectre Xp(k). 10 Cette étape a3 permet également d'obtenir une première estimation O(n) de p la phase de la composante sinusoïdale de fréquence estimée fûK(N/2), à l'échantillon n, par la formule de déroulement de phase suivante: çop û ^K(n) ù çî) (N/2) + 27r n ).fK(N/2)Te, p Te étant la période d'échantillonnage du signal. 15 Dans la quatrième étape a4, on corrige la première estimation de phase (n) de la composante sinusoïdale de fréquence estimée fùK(N/2), afin d'obtenir une seconde estimation de phase, qui est une estimation de phase corrigée -ç5p n) plus proche de la phase réelle çep n) de cette composante sinusoïdale à l'échantillon n de la trame d'indice p. 20 En effet, la première estimation de phase çe(n) donne de bons résultats pour les signaux stationnaires en fréquence, c'est-à-dire pour les signaux dont les fréquences sont constantes pendant toute la durée de l'intervalle d'analyse. Par exemple si xp(n) est un signal sinusoïdal pur défini de la manière suivante: 25 xp(n) = 4 cos(Çep (n)) avec q).pK (n) = nTe + ça() où A0 et (po sont l'amplitude et la phase du signal xp(n) à l'origine, 9pK (n) est sa phase à l'échantillon n, fK sa fréquence et Te la période d'échantillonnage, 2905489 12 alors la première estimation de phase Op (n) rend la phase réellep (n) avec une bonne précision. En revanche, pour des signaux de fréquences variables pendant la durée de l'intervalle d'analyse, la phase estimée Op n) par la méthode 5 décrite ci-dessus présente un biais. Par exemple si xp(n) est un signal défini de la manière suivante: xp (n) = A0 cos(~pp (n)) avec ~p (n) = 2n-(fK + 2K nTe )nTe + çoo où Ao et (po sont l'amplitude et la phase du signal xp(n) à l'origine, 10 çp (n) est sa phase à l'échantillon n, fK sa fréquence, Te la période d'échantillonnage, et fK la variation de la fréquence fK du signal xp(n) sur la durée de l'intervalle d'analyse, correspondant à une variation de phase du deuxième ordre non 15 négligeable, alors la première estimation de phase çp (n) ne rend pas la phase réelle çpp (n) du signal xp(n) à l'échantillon n de manière satisfaisante. En particulier, lorsqu'on modélise les composantes sinusoïdales d'un signal naturel en utilisant cette première estimation de phase, et qu'on soustrait les 20 composantes ainsi modélisées du signal original, on obtient un bruit résiduel important. Le procédé d'estimation de phase selon l'invention utilise une approximation du biais entre la phase réelle p (n) et la première estimation de phase çn (n) , qui est ainsi définie: 2905489 13 B(vf,dfK)=ùarctan avec (Sf = fK (N/2) ù = où B(8f, dfK) est une approximation du biais entre la phase réelle e (n) et la 5 première estimation de phase çb (n) , 8f représente l'écart au milieu de la fenêtre d'analyse entre la fréquence réelle de la composante sinusoïdale du signal xp(n) correspondant au pic fréquentiel d'indice K, et la fréquence discrète correspondante, fK(N/2) est la fréquence réelle de la composante sinusoïdale du signal xp(n) 10 correspondant au pic fréquentiel d'indice K au milieu de la fenêtre d'analyse, K est l'indice fréquentiel du pic correspondant à la composante sinusoïdale considérée dans le spectre de Fourier du signal xp(n), Fe est la fréquence d'échantillonnage du signal, N est la longueur de la fenêtre d'analyse, 15 AfK est la variation de fréquence de la composante sinusoïdale de fréquence fK(N/2) entre le début et la fin de ladite fenêtre d'analyse, n est l'indice d'échantillonnage de la trame, Te est la période d'échantillonnage du signal, et w(n) est la valeur de la fenêtre d'analyse à l'échantillon n. 20 Afin d'obtenir une estimation corrigées (n) de la phase de la composante sinusoïdale de fréquence estimée ,f (N/2) dans la trame d'indice p, lorsque cette phase suit localement une variation quadratique, il faut donc estimer la variation de fréquence de cette composante sur la durée de N 2905489 14 l'intervalle d'analyse. Cette estimation Af. est faite par exemple à partir de l'analyse des partiels construits aux étapes a5 pour les trames précédentes, c'est-à-dire lors du liage des paramètres des pics extraits aux étapes a3 correspondantes. Dans ce cas la variation de fréquence est estimée par la 5 différence entre les fréquences estimées de deux pics correspondants de deux trames successives, divisée par l'intervalle de temps qui les sépare. D'autres méthodes plus sophistiquées d'estimation de la variation de fréquence sont utilisables, telles que celles qui sont présentées dans l'article de Aaron S. Master et Yi-Wen Liu pour la conférence 2003 "International 10 Conference on Acoustics, Speech, and Signal Processing (ICASSP) de l'organisme IEEE, intitulé "Nonstationary Sinusoidal Modeling with Efficient Estimation of Linear Frequency Chirp Parameters". L'approximation du biais définie plus haut est pré-calculée pour différentes valeurs indicées de différence de fréquence 8f entre la fréquence 15 d'une composante sinusoïdale d'un signal numérique et la fréquence discrète correspondante dans le spectre de Fourier de ce signal, ainsi que pour différentes valeurs indicées de variation de fréquence estimée entre le début et la fin de la fenêtre d'analyse du signal numérique. Ces valeurs sont calculées à l'avance et sont donc réutilisables pour 20 toute itération des étapes al à a5 sur différents signaux numériques non stationnaires dès lors que la fenêtre d'analyse et la fréquence d'échantillonnage restent inchangées. Ces valeurs sont organisées dans une grille de valeurs de biais prédéterminées B(8[u],A[v]) pour la fenêtre d'analyse w(n) et la fréquence 25 d'échantillonnage Fe, u et v étant des indices sous forme de nombres entiers. Les valeurs indicées de différences de fréquence ô[u] sont par exemple régulièrement espacées de 1,56 Hz, et les valeurs indicées de variations de fréquence A[v] sont régulièrement espacées de 1 Hz/s (Hertz/ seconde) sur une plage de variation de fréquence allant de 0 à 8000Hz/s. 2905489 15 Afin d'obtenir une estimation corrigée pK (n) de la phase de la composante sinusoïdale de fréquence estimée fK (Nt 2 dans la trame d'indice p, on effectue à l'étape a4 une interpolation linéaire de valeurs de biais prédéterminées de la grille de valeurs. Pour cela, on calcule la différence de 5 fréquence âf entre la fréquence estimée 7K (N/ 2) à l'étape a3 au milieu de la fenêtre d'analyse, et la fréquence discrète correspondante, par la formule: ôf = fK (Ni 2)ùK Fe soit 25f = fK (Ni 2 ù FK où âf représente l'écart entre la fréquence estimée au milieu de la fenêtre 10 d'analyse de la composante sinusoïdale du signal xp(n) correspondant au pic fréquentiel d'indice K, et la fréquence discrète correspondante, fK(N/2) est la fréquence estimée de la composante sinusoïdale du signal xp(n) correspondant au pic fréquentiel d'indice K au milieu de la fenêtre d'analyse, 15 K est l'indice fréquentiel du pic correspondant à la composante sinusoïdale considérée dans le spectre de Fourier du signal xp(n), Fo est la fréquence d'échantillonnage du signal, FK est la fréquence correspondant à l'indice fréquentiel K dans le spectre Xp(k), 20 et N est la longueur de la fenêtre d'analyse. La variation de fréquence estimée 47K entre le début et la fin de la fenêtre d'analyse ayant été précédemment obtenue, on cherche les indices u et v de la grille de valeurs prédéterminées tels que: ≤. (5f [u +1] et Af[v] A ≤. ,àf [1, +1] 25 On estime ensuite le biais entre la première estimation de phase de l'étape a3, et la phase réelle de la composante sinusoïdale considérée du signal xp(n) au milieu de la fenêtre d'analyse, par interpolation linéaire, suivant les équations suivantes: 2905489 16 B = B'+(B"ùB' 4fK -f[v] `[h+t]-4f[v] avec Sf{u],Af[v])+(B(8f[u+1],A et B'= B -B(8f[u],4f[v]) sf[u+l]ù [u] Sf -8f [u] 5 B"=B(:f 8 [u],Af[v+i]) + (B(rf [u+l],Af[v+1]) - B( âf [u] -'5f[u] , Af[, + 1D) sf[u ')f +1]ù~f[u] Il est à noter que dans ce mode de réalisation on effectue une interpolation linéaire de quatre valeurs de la grille de valeurs prédéterminées, mais il est également envisageable de n'utiliser qu'une valeur de la grille de valeurs, ou plus. De plus d'autres types d'interpolation sont utilisables, comme 10 par exemple une interpolation polynômiale. En outre, dans ce mode de réalisation de l'invention, on estime le biais entre la première estimation de phase et la phase réelle de la composante sinusoïdale à l'aide d'une grille de valeurs de l'expression en arctangente donnée précédemment. Cependant d'autres estimations du biais sont 15 possibles. Par exemple, en variante de réalisation, pour de faibles valeurs d'écarts de fréquences entre la fréquence estimée fK (N / 2) et la fréquence discrète FK correspondante, et pour de faibles valeurs de variation de fréquence estimée 7 , on utilise un développement limité de cette expression en arctangente: 20 B= (2W'c5f + Te i , dfK ) W0 où W;, i étant un indice variant de 0 à 2, est le moment d'ordre i de la fenêtre d'analyse, défini par la formule: N -I étant la valeur de la fenêtre d'analyse à l'échantillon n. 2905489 17 Une fois le biais estimé, l'estimation de phase corrigée en N/2 à l'étape a4 est donnée par la formule: "'K çop N/2) =0p(N/2)+B , où 5 N/2 est l'indice de l'échantillon au milieu de la trame d'indice p, -K çop (N/ 2) est la première estimation de phase en N/2 corrigée par le procédé d'estimation de phase selon l'invention, çe (N/ 2) est la première estimation de phaseen N/2 obtenue à l'étape a3, et i3 est l'estimation du biais obtenue ci-dessus, entre la phase réelle de la 10 composante sinusoïdale du signal xp(n) au milieu de la fenêtre d'analyse et la première estimation de phase de cette composante au milieu de la fenêtre d'analyse. L'estimation de phase corrigée (opK (n) à un échantillon n quelconque de la fenêtre d'analyse est également donnée par la formule: 15 (n) Op (n) + B où -0 (n) est l'estimation de phase corrigée par le procédé selon l'invention à l'échantillon n, "K. çop (n) est la première estimation de phase obtenue à l'étape a3 pour 20 l'échantillon n, et B est l'estimation du biais obtenue ci-dessus. Il est à noter que la formule précédente s'obtient par simple déroulement de phase: K K / 2 2 À (N / 2) +A K (n ù / 2)Te )(n ù N / 2)7; 25 Op (N/2)+B+21r(fK(N/2 + (n - N / 2)Te r / 2)Te 2905489 18 ) = (N/2)+27r(7 .N/2)+ K (nûN/2)T2)Te +B fPK (n) _ B où fK. (N/2) est la fréquence estimée de la composante sinusoïdale du signal 5 xp(n) correspondant au pic fréquentiel d'indice K au milieu de la fenêtre d'analyse, fK est la variation de fréquence estimée entre le début et la fin de la fenêtre d'analyse, et Te est la période d'échantillonnage. 10 Les résultats obtenus par le procédé d'estimation de phase selon l'invention sont illustrés dans le tracé TRA de la figure 4. Ce tracé permet de comparer les performances du procédé selon la présente invention avec les résultats obtenus sans correction de phase ou en utilisant des techniques de 15 correction de phase antérieures décrites dans la demande de brevet FR0553891 et qui n'utilisent pour l'analyse spectrale que la Transformée de Fourier Discrète. Les tests pour effectuer ce tracé ont été réalisés à partir d'un grand nombre d'estimations de phase pour différentes fréquences d'un signal 20 synthétique sinusoïdal échantillonné à 16000Hz, auquel on a ajouté un bruit blanc gaussien centré pour obtenir un rapport signal à bruit variant de -10 à 40 dB (décibels). La plage de variations en fréquence utilisée va de 0 à 8000Hz/s. Les erreurs ainsi obtenues entre les estimations de phase 25 éventuellement corrigées et les phases théoriques correspondantes du signal ont été moyennées sur les différents tests correspondant chacun à une estimation de phase sur un signal bruité de fréquence différente pour une valeur de rapport signal à bruit donnée. 2905489 19 Une grille de valeurs prédéterminées de biais a été remplie préalablement à ces tests. Pour cela, des valeurs de biais B(8[u],A[v]) ont été calculées suivant le procédé selon l'invention pour différentes valeurs indicées 8[u] de différences de fréquence et A[v] de variations de fréquences, celles-ci 5 variant de 0 à 8000Hz/s. La résolution utilisée est: e[u +1] ù [u] = 1,5611z Af [v +1] û Af [v] =1,0Hz/s Pour chaque test, correspondant à une fréquence de signal donnée, à une variation de fréquence donnée et à un rapport signal à bruit donné, on a 10 tout d'abord réalisé les étapes suivantes: - on a découpé le signal bruité en trames de longueur N égale à 512 avec un pas d'avancement M égal à 128, on a effectué ensuite une TFD du signal bruité fenêtré par une fenêtre de Hann, soit w(n), et sur lequel on a effectué une opération de "zero- 15 phasing", - on a effectué ensuite trois estimations de fréquence pour trois trames consécutives par une méthode telle que celle du vocodeur de phase, au milieu de chacune des trois trames analysées, ce qui a permis d'obtenir une estimation ÂfK la variation de fréquence du signal bruité, 20 et une estimation "K (N/2) de la fréquence sur la deuxième trame en son milieu, - on a utilisé ensuite la TFD de la deuxième trame, afin d'obtenir une première estimation de phase au milieu de cette trame en prenant l'argument du pic d'énergie maximale de cette TFD, 25 on a calculé l'erreur entre cette première estimation de phase et la phase théorique du signal non bruité au milieu de cette deuxième trame, c'est-à-dire la valeur: 2 2) où 2905489 20 9p (N/2) est la première estimation de phase du signal bruité au milieu de la deuxième trame, et K (N/2) est la valeur théorique de la phase du signal non bruité au milieu de la deuxième trame. 5 L'erreur ainsi calculée a permis d'établir la courbe CSC d'erreurs quadratiques moyennes d'estimations de phase en fonction du rapport signal à bruit lorsqu'on n'utilise aucune technique de correction de la première estimation de phase. Consécutivement à ces premières étapes, on a ensuite appliqué, pour chaque 10 test, différentes techniques de correction de la première estimation de phase: - on a corrigé la phase estimée précédemment suivant une première variante du procédé de correction de phase décrit dans la demande FR0553891, par la formule: (N/2) =0pK(N/2)-KT,2((W"2 ùM'Vl) +N2)' Wo 4 15 où N est la longueur de la fenêtre d'analyse, Te est la période d'échantillonnage, 41K est la variation de fréquence précédemment estimée du signal bruité, 20 0,K(N/2) est la première estimation phase du signal bruité sur la deuxième trame en son milieu, (N/ 2) est la phase estimée corrigée du signal sur la deuxième trame en son milieu, et Wi, i étant un indice variant de 0 à 2, est le moment d'ordre i de la 25 fenêtre d'analyse, défini par la formule: o 2905489 21 w(n) est la valeur de la fenêtre d'analyse à l'échantillon n. et on a calculé l'erreur entre la phase estimée corrigée précédemment et la phase théorique du signal non bruité au milieu de la deuxième trame, c'est-à-dire la valeur: 5 (N / 2) - (N / 2) L'erreur ainsi calculée a permis d'établir la courbe CL d'erreurs quadratiques moyennes d'estimations de phase corrigées suivant la première variante du procédé de correction de phase de la demande FR0553891, en fonction du rapport signal à bruit. 10 On a ensuite corrigé la première estimation de phase suivant une seconde variante du procédé de correction de phase décrit dans la demande FR0553891, par la formule: i=o avec Çà,K (N / 2) = (N/ 2) aiAfK i - n-A:i K.Te2 N 4 15 d=4 ao =9,3599.10- al =1,0145.104 a2 =1,9666.10-9 a3 = -1,0262.10-12 20 a4 = 4,4716.10-17 , où N est la longueur de la fenêtre d'analyse, Te est la période d'échantillonnage, AfK est la variation de fréquence précédemment estimée du signal 25 bruité, vp (N/2) est la première estimation phase du signal bruité sur la deuxième trame en son milieu, 2905489 22 et û3p N/2) est la phase estimée corrigée du signal sur la deuxième trame en son milieu suivant cette deuxième variante du procédé de correction décrit dans la demande FR0553891, et on a calculé l'erreur entre cette nouvelle phase estimée corrigée 5 'e, (N/2) et la phase théorique du signal non bruité au milieu de la deuxième trame, c'est-à-dire la valeur: (N 2)-çpp (N/2) L'erreur ainsi calculée a permis d'établir la courbe CP d'erreurs quadratiques moyennes d'estimations de phase corrigées suivant la 10 deuxième variante du procédé de correction de phase de la demande FR0553891, en fonction du rapport signal à bruit. Enfin on a corrigé la première estimation de phase suivant le mode de réalisation du procédé selon l'invention par la formule: ço` (N/2) = Op' (N / 2) +1-3 avec B = B'+(B"ûB') dfK û 4f [v] B'= B(ff'[u],4f[v ) B(Sf[u+1],of[v])-B(sf[u],.f[v])) -8f[u] et 8f -8f[u] )Sf[u+1]-âf[u] u et v étant les indices des valeurs de biais B(6[u],[v])de la grille de 20 valeurs prédéterminées, tels que: 8f [u] bf s 8f [u + 1] et f [v] .≤ LfK REVENDICATIONS 1. Procédé d'estimation de phase d'une composante sinusoïdale d'un signal numérique (x(n)) non stationnaire, comportant les étapes de: segmentation (a1) dudit signal en trames temporelles, - Transformée de Fourier Discrète (a2) d'une desdites trames, pondérée par une fenêtre d'analyse (w(n)) afin d'obtenir une valeur d'un pic du spectre en fréquence dudit signal correspondant à ladite composante sinusoïdale, - estimation de fréquence (a3) sur ladite trame de ladite composante sinusoïdale, - première estimation de phase (a3, (n)) de ladite composante sinusoïdale en utilisant l'argument de ladite valeur de pic ainsi obtenue, caractérisé en ce qu'il comprend les étapes supplémentaires de: - estimation (a4) d'une valeur de biais en fonction de la différence de fréquence entre la fréquence précédemment obtenue à l'étape d'estimation et la fréquence discrète dudit pic, et en fonction d'une variation de fréquence estimée de ladite composante sinusoïdale ( AiK ) entre le début et la fin de ladite fenêtre d'analyse (w(n)), - correction (a4) de ladite première estimation de phase par ladite valeur de biais. 2. Procédé d'estimation de phase d'une composante sinusoïdale d'un signal numérique (x(n)) non stationnaire selon la revendication 1, caractérisé en ce que l'étape d'estimation d'une valeur de biais utilise une grille de valeurs prédéterminées( B(8[u],A[v] comprenant, pour des valeurs de différences de fréquences et de variations de fréquence prédéterminées, des valeurs de biais 2905489 26 B 8 ], ] = ùaretan où B(8[u], A[v]) est une valeur de biais de la grille de valeurs prédéterminées, 8[u] est une valeur de différence de fréquences prédéterminée, 5 A[v] est une valeur de variation de fréquence prédéterminée, n est l'indice d'échantillonnage de la trame, N est la longueur de la trame, Te est la période d'échantillonnage du signal, et w(n) est la valeur de la fenêtre d'analyse à l'échantillon n. 1 0 3. Procédé d'estimation de phase d'une composante sinusoïdale d'un signal numérique (x(n)) non stationnaire selon la revendication 2, caractérisé en ce que dans l'étape d'estimation (a4) d'une valeur de biais, ladite valeur de biais est obtenue par interpolation de valeurs prédéterminées de ladite grille. 4. Procédé de modélisation harmonique d'un signal numérique (x(n)) non stationnaire, utilisant le procédé d'estimation de phase selon l'une quelconque des revendications 1 à 3. 20 5. Dispositif d'estimation de phase d'une composante sinusoïdale d'un signal numérique (x(n)) non stationnaire, caractérisé en ce qu'il comporte des moyens de mise en oeuvre du procédé selon l'une quelconque des revendications à 3. 15 2905489 27 6. Programme d'ordinateur caractérisé en ce qu'il comporte des instructions de mise en oeuvre du procédé selon l'une quelconque des revendications 1 à 3, lorsque ledit programme est exécuté dans un système informatique.