PROCEDE ET DISPOSITIF PERMETTANT LE SUIVI DE DOPPLER POUR

MODEM A LARGE BANDE

L'invention concerne un procédé et un dispositif permettant le suivi Doppler, pour modem à large bande.

Plus généralement, elle concerne aussi un système d'estimation et de suivi du décalage en fréquence pour un système de transmission dont la largeur de bande n'est pas « négligeable » devant la fréquence centrale.

Le problème posé est la transmission de données de type quelconque entre deux correspondants immergés dans l'eau et à un débit de plusieurs centaines de bits par seconde tout en utilisant une bande passante aussi faible que possible et une puissance de calcul réduite au minimum, en particulier pour des équipements autonomes fonctionnant sur batteries.

Le problème est aussi celui de la transmission de données entre deux correspondants immergés dans un milieu donné.

Les ondes électromagnétiques classiquement utilisées (de quelques centaines de kHz à quelques GHz se propageant extrêmement mal dans l'eau, le modem (modulateur-démodulateur) utilise donc des ondes acoustiques pour véhiculer les messages transmis.

Pour rappel, le milieu sous-marin a un comportement extrêmement complexe pour les ondes acoustiques. Ses inhomogénéités (salinité, température, courants, ..) font que le chemin suivi par les ondes n'est une ligne droite qu'à très courte distance et que ce trajet varie au cours du temps.

Les ondes qui arrivent au récepteur peuvent suivre plusieurs trajets avec des délais très différents. Ces trajets peuvent être soit « directs » (plusieurs parcours différents convergeant vers le récepteur), soit « indirects » (réflexions sur la surface de l'eau, éventuellement agitée de vagues, ainsi que sur le fond qui peut avoir toutes sortes de types de reliefs).

Dans la pratique, cela se traduit par l'arrivée au récepteur de nombreuses répliques du signal émis, chacune ayant à un instant donné, une amplitude, une phase, un retard et un décalage fréquentiel propre. A ceci, il faut ajouter, éventuellement, un décalage fréquentiel du à la vitesse radiale relative des correspondants (effet Doppler) ainsi qu'un décalage

fréquentiel pseudo périodique dû à la houle dans le cas où l'un au moins des correspondants est à la surface.

Dans certains cas, la fiabilité de la transmission doit être quasi parfaite.

S'il s'agit d'un texte, on peut tolérer quelques erreurs isolées, car elles seront automatiquement corrigées par le lecteur. Lorsqu'il s'agit d'un fichier binaire (par exemple une image JPEG), le taux d'erreurs exigé est voisin de 0, sous peine de très forte dégradation de l'image.

De plus, le système est « à large bande ». Dans les applications pratiques, la largeur de bande du signal émis n'est pas négligeable comparée à la fréquence centrale d'émission. Le décalage Doppler étant proportionnel à la fréquence, les fréquences basses ont un décalage en fréquence plus faible que les fréquences hautes.

Qui plus est, la vitesse relative de l'émetteur et du récepteur n'étant pas toujours négligeable par rapport à la vitesse du son dans l'eau (environ 1500 m/s) on observe un effet de dilatation ou de contraction du temps selon le signe de la vitesse relative.

Par exemple, un message émis pendant une durée de 10s sera reçu pendant 10.13 secondes si la vitesse d'éloignement est de +10 m/s (36 km/h ou environ 20 nœuds) ce qui pose des problèmes de suivi de synchronisation.

A la connaissance du demandeur, il existe des systèmes de type « à étalement de spectre » dans lesquels le débit utile est très faible comparé à la bande passante utilisée et où la discrétion est de mise, ce qui interdit l'utilisation de porteuses fixes, facilement repérables. La demande de brevet WO 2005/055442 divulgue un système « à large bande » qui estime une seule fois et pour toute la durée du message une valeur de décalage Doppler sans effectuer de suivi.

L'objet de la présente invention concerne notamment un procédé permettant d'estimer en continu le décalage en fréquence (Doppler) pour un modem fonctionnant en milieu sous-marin.

L'invention concerne notamment un système d'estimation et de suivi du décalage en fréquence pour un système de transmission dont la largeur de bande

n'est pas négligeable devant la fréquence centrale, le système émettant une forme d'onde incluant des porteuses fixes de niveau sensiblement supérieur à celui des porteuses véhiculant l'information, les porteuses ayant des fréquences Fo ... FN-1 et le signal utile porteur de l'information caractérisé en ce qu'il comporte au moins les éléments suivants :

• N bancs de filtres passe-bande, chaque banc comportant P filtres, adaptés à couvrir l'ensemble des valeurs que peuvent prendre les fréquences Fo ... FN-1 après décalage Doppler,

• La sortie d'un filtre est reliée à un détecteur quadratique d'enveloppe, • P sommateurs adaptés à additionner les sorties correspondantes des N bancs de filtres,

• Un système de détection de la position des deux plus grands maxima parmi les P sommes et de validation de la position du plus grand maxima,

• Un système de mémorisation des positions des deux plus grands maximas et de leurs « poids » et de prédiction de la valeur courante,

• Un système de validation/remplacement de la position trouvée,

• un système convertissant les positions en Doppler estimé et filtrant celui-ci par un filtre passe-bas.

L'invention concerne aussi un procédé permettant d'estimer et de suivre l'effet Doppler pour un signal composé d'une forme d'onde incluant des porteuses fixes de niveau sensiblement supérieur à celui des porteuses véhiculant l'information, les porteuses ayant des fréquences Fo ... FN-1 et du signal utile porteur de l'information dans un système de transmission dont la largeur de bande n'est pas négligeable devant la fréquence centrale caractérisé en ce qu'il comporte au moins les étapes suivantes :

1 . Transmettre le signal composite comportant le signal utile porteur de l'information et le signal de référence à une étape de filtrage utilisant les fréquences des porteuses,

2. Déterminer pour chaque filtre, l'enveloppe du signal, e _n , _p , 3. Sommer les différentes sorties e _n , _p ,

4. Déterminer les valeurs des deux plus grands maximas, M et m, pour plusieurs instants, t,

5. Déterminer la valeur probable de la position à l'instant t, en considérant les maximas trouvés qui définissent de façon unique la valeur de décalage Doppler associée,

6. Effectuer un filtrage passe-bas sur les valeurs de décalage Doppler pour obtenir une valeur lissée.

Le procédé comporte par exemple une étape où l'on interpole la valeur lissée à des instants intermédiaires au moyen de filtres passe-bas dérivés du filtre défini précédemment.

L'étape de filtrage 1 comporte, par exemple, au moins les étapes suivantes :

• On effectue le fenêtrage du signal d'entrée par une fenêtre de Hanning (simple cosinus sur-élevé) sur une durée de 160 ms, avec recouvrement de 80 ms entre deux fenêtres successives,

• aux 1 .024 points obtenus, on ajoute 3.072 points nuls, puis : « on effectue la transformée de Fourier rapide FFT 4096 points sur ce signal dont la sortie est notée Z ₀ , Z ₁ ... Z ₄₀₉₅

• pour une fréquence centrale F _n , _p souhaitée, on calcule le rang n du "bin" Z _n de la transformée de Fourier FFT qui est juste inférieur,

• on calcule le nombre réel x compris entre 0 et 1 qui vaut (F _n , _p - 1 .5625 n) / 1.5625 et l'on sélectionne deux nombres complexes zθ(x) et z1 (x) fonctions de x,

• on détermine la sortie du filtre par zθ(x) Zn + z1 (x) Zn+1.

Le procédé et le système selon l'invention présentent notamment les avantages suivants. Ils permettent d'effectuer des communications sous-marines de longue durée, à Doppler variable, entre un émetteur et un récepteur en déplacement sur des trajectoires quelconques.

• La figure 1 , dans un espace temps-fréquence, l'effet du décalage sur les porteuses,

• La figure 2, un exemple de structure des modules et des étapes mises en œuvre par le système selon l'invention, « La figure 3, une étape de filtrage du signal par fenêtrage du signal d'entrée,

• La figure 4, un exemple de sorties des détecteurs,

• La figure 5, un exemple des sorties des sommateurs,

• La figure 6, un exemple de régression parabolique,

• Les figures 7 à 9, les résultats du filtrage final, et « La figure 10, un exemple de réponse en fréquence obtenue.

Afin de mieux faire comprendre le principe mis en œuvre dans l'invention, l'exemple qui suit est donné pour un modem fonctionnant en milieu sous- marin, caractérisé par un étalement temporel important et rapidement variable et une largeur de bande non négligeable devant la bande passante qui rend le décalage en fréquence du à l'effet Doppler variable dans la bande passante même du signal (effet « bande large »).

La forme d'onde émise par le modem en question inclut des porteuses fixes (non modulées) de niveau sensiblement supérieur à celui des porteuses véhiculant l'information. C'est à partir de ces porteuses que se fait l'estimation du décalage Doppler moyen du signal, décalage susceptible de varier rapidement.

Les porteuses de référence sont en nombre égal à N et ont des fréquences notées F ₀ , F ₁ , ... F _N . _! .

Le décalage Doppler maximum est donné par un nombre sans dimension, noté δ, obtenu par l'équation : _{δ =} V _r (R) - V _r (E) C - V r (R)

C : vitesse du son dans l'eau v _r (E) : vitesse radiale de l'émetteur v _r (R) : vitesse radiale du récepteur Référentiel : l'eau, supposée immobile

En présence d'un tel décalage Doppler, chaque porteuse de fréquence F, voit sa fréquence modifiée en (1 + δ) F ₁ .

La figure 1 décrit dans un espace temps-fréquence, l'effet du décalage Doppler sur les porteuses, pour δ = -δ _max , d = 0, et d = +δ _max , δ _max étant un élément de conception du système.

Les fréquences choisies peuvent être quelconques, mais, afin que l'estimation de δ puisse se faire sans ambiguïté, il faut que les petits rectangles 1 de la figure 1 aient une hauteur positive, ce qui se traduit par la condition simple suivante :

(1 - δ _max ) F _n > (1 + δ _max ) F _n _ ! ou

Fn - Fn - I > δ _max ( F _n + F _n _ i ) n = 1 ... N-I

La figure 2 représente les différentes étapes mises en œuvre dans l'invention et les éléments constituant le dispositif.

Ce système nécessite la présence, pendant toute la durée du message et même un peu plus, d'un ensemble de N porteuses fixes non modulées en sus du signal utile. Le signal complexe (porteuses de référence et signal utile) est reçu sur une chaîne de traitement comportant, dans cet ordre :

- un ensemble, 100, de N bancs de P filtres passe-bande couvrant (légèrement plus que) au moins l'ensemble des valeurs que peuvent prendre les fréquences F ₀ ... F _N . _! après décalage Doppler. Chaque filtre passe-bandes est suivi d'un détecteur quadratique d'enveloppe non représenté distinctement des filtres sur la figure pour des raisons de simplification.

- un ensemble 200, de P sommateurs additionnant les sorties des détecteurs quadratiques d'enveloppe de même rang pour chaque filtre, - un système 300, de détection de la position des deux plus grands maxima parmi les P sommes et de validation de leurs positions en sortie de la série de sommateurs,

- un système 400 de mémorisation des positions et de leurs « poids » et de prédiction de la valeur courante, un système 500 de validation / remplacement de la position trouvée,

- un système 600 convertissant les positions en Doppler estimé et filtrant celui-ci par filtre bas qui donne une estimation du décalage Doppler

(avec un retard fixe).

Le mécanisme d'estimation et de suivi peut aussi inclure un mécanisme d'interpolation permettant de disposer plus fréquemment d'échantillons de décalage Doppler.

Le détail des différents éléments faisant partie du système selon l'invention est donné ci-après.

Bancs de filtres 100

Un banc de filtres est adapté notamment à couvrir les fréquences comprises entre (1 - δ _max ) F _n et (1 + δ _max ) F _n où F _n est la fréquence après décalage Doppler pour le banc de filtre de rang n.

Chaque banc de filtres comporte, par exemple, le même nombre P de filtres passe-bande de fréquences centrales régulièrement espacées.

Le p ^eme filtre du n ^eme banc sera centré sur une fréquence égale à :

2 p - P + 1

F ^r n, p — — F ^r n H + δ ',max n = 0 ... N-l p = 0 ... P-I

F _n : la fréquence après décalage Doppler pour le banc de filtre de rang n. Une manière simple de réaliser le banc de filtres est décrite ci-dessous.

A titre d'exemple, on suppose que les porteuses sont espacées d'un multiple de 25 Hz et que la fréquence d'échantillonnage est de 6.400 Hz.

On suppose que, lors de la phase de démodulation, un nouveau symbole apparaît toutes les 80 ms et que l'on a alors besoin d'une nouvelle valeur du décalage en fréquence.

En supposant que la démodulation de phase différentielle peut avoir jusqu'à 4 états de phase, soit au minimum 90 °. Si l'erreur d'estimation du décalage en fréquence vaut dF, la rotation de phase supplémentaire qu'il introduit est de 360° x 80 ms x dF.

Pour qu'elle ne dépasse pas le quart de la rotation minimum (règle empirique issue de l'expérience), l'erreur d'estimation du décalage en fréquence dF doit valoir au maximum df _MA x = 0.78125 Hz (en valeur absolue), ce qui signifie que la transformée de Fourier rapide, FFT, à utiliser doit pouvoir séparer deux fréquences séparées de 2 df _MA x = 1.5625 Hz. Son nombre de points devrait alors être d'au moins 6.400/1.5625 = 4.096.

Pour cela, on procède de la manière suivante : « On effectue le fenêtrage (pondération) du signal d'entrée par une fenêtre de

Hanning (simple cosinus sur-élevé) sur une durée de 160 ms, avec recouvrement de 80 ms entre deux fenêtres successives. Cette étape est illustrée à la figure 3.

• aux 0.16 x 6.400 = 1 .024 points obtenus, on ajoute 4.096 - 1 .024 = 3.072 points nuls, puis :

1. on effectue la transformée de Fourier rapide FFT 4096 points sur ce signal dont la sortie est notée Z ₀ , Z ₁ ... Z ₄₀₉₅ ,

2. pour une fréquence centrale F _n , _p souhaitée, on calcule le rang n du "bin" Z _n de la transformée de Fourier FFT qui est juste inférieur, c'est à dire celui pour lequel F _n , est comprise entre n fois 1 .5625 Hz et (n + 1 ) fois

1.5625 Hz,

3. on calcule le nombre réel x compris entre 0 et 1 qui vaut (F _n , _p - 1 .5625 n) / 1.5625 et l'on sélectionne deux nombres complexes zθ(x) et z1 (x) fonctions de x (et pré-calculés par commodité), 4. la sortie du filtre est simplement donnée par zθ(x) Zn + z1 (x) Zn+1.

Les coefficients zO et z1 sont optimisés afin que le résultat soit aussi proche que possible de celui que donnerait un filtre réellement centré sur F _n , _p .

L'un des deux coefficients zO ou z1 peut être choisi REEL (et l'autre corrigé en conséquence) puisque la phase du signal de sortie est sans importance. Chaque filtre est suivi d'un détecteur quadratique d'enveloppe donnant le carré du module de sa sortie (le traitement est généralement fait sur un signal complexe, en bande de base, i.e. centré sur la fréquence O).

La sortie du p ^eme détecteur du n ^eme banc sera notée e _n,p .

La figure 4 donne un exemple de sorties des détecteurs quadratique d'enveloppe, pour P = 100, et N = 7 fréquences, avec un signal d'entrée décalé vers les hautes fréquences. Sommateurs 200

La sortie du p ^eme sommateur est donnée par :

N-I

Sp = ∑e _n , _P

B=O p = 0....P - I avec e _n , _p .la sortie du p ^eme détecteur du n ^eme banc de filtre.

La figure 5 représente un exemple des sorties des sommateurs pour le cas représenté sur la figure 4, où l'on voit bien l'effet de lissage apporté par la sommation (du fait de la propagation, le gain peut varier beaucoup au voisinage de chacune des N porteuses). Détection des deux plus grands maxima parmi P 300

Le signal issu des sommateurs 200 est transmis à un système de détection des maximas 300 ayant notamment pour fonction de rechercher Deux maxima.

En effet, s'il existe deux maxima d'amplitudes voisines, ils sont en général de faible amplitude et ne correspondent à rien de précis : il y a probablement uniquement du bruit en entrée du récepteur.

Pour la recherche de ces deux maximas, M (le plus grand) et m (le second plus grand), le procédé peut utiliser l'algorithme suivant, écrit en langage C simplifié.

/ ^* Initialisation 7 M = O position M = 0 m = 0 position m = 0

/ ^* Boucle pour rechercher un maximum aux positions 1 ... P-2 7 pour( i = 1 à P - 2 ) { si ( S(i) >= S(i - 1 ) et S(i) >= S(i + 1 ) ) { / ^* S(i) est un maximum local 7 si ( S(i) > M ) {

/ ^* S(i) est le nouveau plus grand maximum 7 m = M position m = position M M = S(i) position M = i } sinon { si ( S(i) > m ) {

/ ^* S(i) est le nouveau second maximum 7 m = S(i) position m = i

} }

} }

Si position M est non nul (0 est une valeur illégale) et si m est inférieur à une proportion donnée de M (par exemple 50 %), M et sa position (position M) peuvent être considérés comme valables. On a trouvé 2 maximas.

Sinon ils sont ignorés par la suite (on se contente de mettre M à zéro).

La valeur de la position des deux maximas est ensuite affinée par exemple par régression quadratique, c'est-à-dire que les amplitudes des sommes

voisines de la somme maximale M sont approximées par une parabole utilisant 3 ou 4 voisins selon le nombre de voisins disponibles, valeur désignée par la suite valeur corrigée.

Pour l'élaboration de l'algorithme de recherche, on pourra utiliser tout autre langage que le langage C.

La figure 6 donne un exemple de régression parabolique, toujours sur les mêmes données des figures 4 et 5.

La position du maximum de la parabole est la valeur modifiée de position M qui, d'entier, devient un nombre réel.

Pour position M = 1 , la correction est égale à :

11 S ₀ - 3 S ₁ - 7 S ₂ - S ₃ 10 ( S ₀ - S ₁ - S _{2 +} S ₃ ) Pour position M = 2....P-3, la correction est égale à :

7(2S _p _ _{2 +} S _p _ _{1 +} S _p+1 - 2S _p+2 ) lθ(2S _p _ ₂ - S _p _ ₁ - 2S _p - S _p+1 + 2S _p+2 )

Pour position_M=P-2 la correction est égale à :

En sortie du détecteur de maximas on dispose de deux valeurs de maximas et de leur position.

Prédicteur de position 400

Le système de prédiction / remplacement, 400, utilise la valeur de la position corrigée à l'instant t, et aux instants t - 1 , t - 2, ..., t - T (t est pris entier pour la commodité des calculs et au moins égal à 4). La valeur de la position corrigée correspond à la valeur de la position déduite du maxima et affinée par régression quadratique.

Sa fonction est notamment d'évaluer, à partir du passé, la valeur probable de la position à l'instant t, en considérant d'abord les maximas de plus grande amplitude (M _M , M, _-2 , ..., M,. _τ ), car les plus fiables. Il est constitué, notamment d'une ligne à retard contenant les T+1 dernières paires M/position_M. Sachant que, sur une durée limitée, la variation du Doppler peut être assimilée localement à une parabole, on prédit la valeur de la position à l'instant t au moyen d'une régression quadratique pondérée.

Pour cela, on calcule les coefficients a, b, c d'une parabole minimisant l'erreur quadratique pondérée suivante (où la variable pos, désigne la valeur de position M après correction) :

E = ∑ Mi((ai ² + bi + c) - posi) ² ι=t-T

La valeur prédite est simplement

pos _t = a t + b t + c Son calcul est simple et donné par la série d'équations suivantes

U _k = ∑M _ι (t -i) ^k ,k = 0 4 ι=t-T

W _k = ∑M _ιP o _Sι (t - i) ^k ,k = 0 2 ι=t-T

_ W ₀ (U ₂ U ₄ - Ul) + W ₁ (U ₂ U ₃ - U ₁ U ₄ ) + W ₂ (U ₁ U ₃ - U ₂ ² )) pos _t = U _{) {U ₂ U ₄ - U ₃ ² ) + U ₁ (U ₂ U ₃ - U ₁ U ₄ ) + U ₂ {U ₁ U ₃ - UI)

Choix de la position à utiliser, 500

La valeur prédite pour la position est transmise à un dispositif 500 de choix final de la position à utiliser. Ce dernier reçoit aussi l'information issue de l'étape de validation/affinage de la valeur de position.

Si la position corrigée est non valide, le calcul est fait avec M, = 0, et l'on remplace alors la valeur prédite « pos, » par sa valeur estimée au moyen de la formule donnée ci-dessus et M, par la moyenne des valeurs précédentes.

Si la position corrigée est valide, on calcule sa valeur estimée avec la valeur mesurée de M,. Si la différence entre la position estimée et la position corrigée est suffisamment faible (en valeur absolue), on garde la valeur de pos, sinon, on fait comme il a été explicité précédemment.

Ce mécanisme permet d'éliminer les valeurs aberrantes.

L'initialisation se fait avec T + 1 valeurs M, égales et petites et des valeurs pos, correspondant à un décalage en fréquence nul, soit (P - 1 ) / 2.

L'approximation de la variation du Doppler par une sinusoïde serait peut- être plus appropriée au canal acoustique sous-marin (à cause du caractère périodique de la houle. Les calculs sont cependant plus complexes.

La conversion de position en "Doppler" peut se faire après la correction de la position, simplement en remplaçant pos, par :

Les calculs seront les mêmes, simplement en remplaçant pos, par δ,.

2 pos , - P + 1 1

O ₁ = δ _π Filtrage final, 600

Ce module reçoit la valeur issue de la conversion de position en Doppler. Il comporte par exemple plusieurs filtre passe-bas et un commutateur. La dernière étape du traitement consiste à effectuer un filtrage passe- bas conséquent sur les valeurs des décalages Doppler pour obtenir une valeur lissée finale δ aussi fidèle que possible tant que la fréquence instantanée de δ reste comprise entre deux fréquences -F ₀ et + F ₀ .

On utilisera par exemple un filtre passe-bas à phase linéaire avec un nombre impair de coefficients, à savoir :

avec ι=2K

!=0 et

(2)h _ι = h _2K _ _ι i = 0....K -l

La condition (1 ) garantit qu'un décalage constant en entrée sera inchangé en sortie.

La condition (2) exprime la phase linéaire, qui garantit un retard constant quelle que soit la vitesse de variation instantanée de δ. Ce retard est égal à K fois la période d'échantillonnage. Il faut lui ajouter le retard dû aux traitements précédents, notamment l'étape de filtrage pour laquelle tous les N x P filtres doivent avoir le même retard, ce qui est le cas avec l'exemple de mise en œuvre proposé.

Pour le design exact du filtre, qui doit être aussi "plat" que possible entre - F ₀ et + F ₀ , on peut utiliser la méthode suivante : 1 - se donner un filtre idéal dont la réponse en fréquence est donnée par :

h(f) = 1 + α -*-

F ₀

- F ₀ <= f <= +F ₀ 0 ailleurs

un exemple est donné à la figure 7 pour F ₀ = 1 Hz et α = 0.2.

2 - calculer sa réponse impulsionnelle (modifiée pour valoir 1 pour t = 0) θ = 2 _π F ₀ t 2 α θ cos( θ) + ( θ ² (1 + α) - 2 α) sin( θ) h(t) = 3

(3 + α) θ ³

(avec : h(0) = 1) donnée par : un exemple de réponse impulsionnelle est donné à la figure 8.

3 - calculer sa réponse impulsionnelle échantillonnée à des intervalles égaux à T, et la pondérer par une fenêtre de Hamming, par exemple : h _{1 =} h ((i- K) T ) 0.54 + 0.46 cos π i - K

K + 0.5 i = 0 ... 2 K

Pour T = 80 ms, et K = 16, par exemple, on a une réponse représentée à la figure 9. 4 - ajuster le coefficient α pour avoir une réponse en fréquence aussi plate que possible dans la bande,

5 - normaliser de telle sorte que la somme des coefficients vaille 1 .

On a représenté en figure 10 un exemple de réponse en fréquence obtenue (la valeur de α est volontairement un peu trop élevée).

Dans le cas où l'on désirerait avoir des valeurs du Doppler à une cadence plus rapide que toutes les 80 ms, on doit utiliser une technique d'interpolation pour obtenir les valeurs intermédiaires.

La technique bien connue de l'homme de l'art et dénommée "filtres polyphasés" permet d'obtenir les échantillons intermédiaires.

Dans un exemple de mise en œuvre, les coefficients h', d'un filtre donnant des échantillons décalés en temps sont calculables simplement par la formule générique :

π p - θ ~γ p + i - θ p - i - θ

suivante, où θ dépend du retard souhaité.

La valeur θ = 0 correspond au filtre d'origine.