FR2724028A1 - Procede d'estimation aveugle de retards differentiels entre deux signaux - Google Patents
Procede d'estimation aveugle de retards differentiels entre deux signaux Download PDFInfo
- Publication number
- FR2724028A1 FR2724028A1 FR9410325A FR9410325A FR2724028A1 FR 2724028 A1 FR2724028 A1 FR 2724028A1 FR 9410325 A FR9410325 A FR 9410325A FR 9410325 A FR9410325 A FR 9410325A FR 2724028 A1 FR2724028 A1 FR 2724028A1
- Authority
- FR
- France
- Prior art keywords
- sep
- signals
- alpha
- vectors
- delays
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
- G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
- G01S3/802—Systems for determining direction or deviation from predetermined direction
- G01S3/808—Systems for determining direction or deviation from predetermined direction using transducers spaced apart and measuring phase or time difference between signals therefrom, i.e. path-difference systems
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
- G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
- G01S3/86—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves with means for eliminating undesired waves, e.g. disturbing noises
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
L'invention concerne les procédés de traitement en aveugle des signaux reçus de deux sources distinctes inconnues, dans une même bande de fréquence. Elle consiste après avoir effectué une transformée de Fourier sur ces signaux, à déterminer leurs multispectres (103), puis les coefficients d'amortissement (104, 105) de ces signaux et enfin les retards (106 - 110) entre ces signaux. Elle s'affranchit du bruit ambiant gaussien, et permet d'éliminer un brouilleur localisé. Les applications vont aux radars et sonars.
Description
La présente invention se rapporte aux procédés qui permettent de déterminer de manière aveugle, c'est-àtire sans rien connaître d'autre, le retard différentiel entre deux signaux provenant de deux sources distinctes.
L'estimation d'un temps de retard différentiel est un vieux problème en traitement du signal. Une collection d'articles de synthèse peut être trouvée dans le numéro spécial de IEEE Trans. ASSP, vol.29, june 1981. Dans ces travaux ainsi que dans les articles de Y.T. CHAN, J.M.
RILEY, and J.B. PLANT "A parameter estimation approach to time-delay estimation and signal detection" IEEE Trans ASSP, 28(1): 8-16, February 1980 l4], de C. NIKIAS and R. PAN "Time delay estimation in unknown gaussian spatially correlated noise" IEEE Trans. ASSP, 36(11): 1706-1714,
November 1988 [8], et de J. TUGNAIT "On time delay estimation with unknown spatially correlated gaussian noise using fourth-order cumulants"
IEEE Trans. Sig. Proc, 39(6): 1258-1267,1991 [9], il n'est question que d'un seul retard. Autrement dit, on suppose que s2(t) = O. Les seuls travaux ayant abordé le problème de l'estimation simultanée aveugle (c. à d. où l'on ne connaît pas les signaux sj(t) ni leurs spectres) de deux retards sont ceux de
J.F. CARDOSO "Source separation using higher order moments "In Proc.
November 1988 [8], et de J. TUGNAIT "On time delay estimation with unknown spatially correlated gaussian noise using fourth-order cumulants"
IEEE Trans. Sig. Proc, 39(6): 1258-1267,1991 [9], il n'est question que d'un seul retard. Autrement dit, on suppose que s2(t) = O. Les seuls travaux ayant abordé le problème de l'estimation simultanée aveugle (c. à d. où l'on ne connaît pas les signaux sj(t) ni leurs spectres) de deux retards sont ceux de
J.F. CARDOSO "Source separation using higher order moments "In Proc.
ICASSP Glasgow, pages 2109-2112, 1989 [2], et de P. COMON "Separation of sources using high-order cumulants" In SPIE Conference on Advanced
Algorithms and Architectures for Signal Processing, pages 170-181, San
Diego, Califomia, August 8-10 1989. vol. Real-time signal processing XII [6], ainsi que ceux listés dans l'article de P. COMON "Independant Component
Analysis, a new concept ?" Signal Processing, 36(3), 1994 [7]. Les approches proposées sont basées sur l'utilisation en cascade des moments d'ordre 2 et des cumulants d'ordre 4 dans le domaine spectral. Ceci présente l'inconvénient d'être sensible à la présence de bruit Gaussien lorsque le rapport signal à bruit est faible. Par ailleurs, dans l'article de J.F.
Algorithms and Architectures for Signal Processing, pages 170-181, San
Diego, Califomia, August 8-10 1989. vol. Real-time signal processing XII [6], ainsi que ceux listés dans l'article de P. COMON "Independant Component
Analysis, a new concept ?" Signal Processing, 36(3), 1994 [7]. Les approches proposées sont basées sur l'utilisation en cascade des moments d'ordre 2 et des cumulants d'ordre 4 dans le domaine spectral. Ceci présente l'inconvénient d'être sensible à la présence de bruit Gaussien lorsque le rapport signal à bruit est faible. Par ailleurs, dans l'article de J.F.
CARDOSO "Iterative techniques for blind source separation using only fourth order cumulants" In Proc. EUSIPCO, Brussels, Belgium, 1992 [3], une approche basée uniquement sur les cumulants d'ordre 4 est proposée, mais cette dernière est itérative, et sa convergence n'est pas assurée.
La présente invention propose une solution analytique, nécessitant donc un nombre fini d'opérations. Cette invention s'affranchit donc des problèmes de convergence. La différence essentielle par rapport à la méthode proposée par E. CHAUMETTE, P. COMON, et D. MULLER "An
ICA-based technique for radiating sources estimation; application to airport surveillance" IEE Proceedinds - Part F, 140(6):39501, December 1993.
ICA-based technique for radiating sources estimation; application to airport surveillance" IEE Proceedinds - Part F, 140(6):39501, December 1993.
Special issue on Applications of High-Order Statistics [5], conceme une réduction considérable de la complexité, et la possibilité de se passer des moments d'ordre 2, ce qui confère au présent procédé une insensibilité potentielle totale aux bruits Gaussiens pour des temps d'intégration suffisamment longs.
Le probléme peut se poser de la maniére analytique suivante: On suppose que l'on reçoit sur deux capteurs les signaux à temps discret x1(t) et x2(t) suivants:
xj(t) = s,(t) + B2 s2(t-#2) + v1(t), (1)
x2(t) = s2(t) + B1 s1(t-#1) + v2(t), (2)
où si (t) sont des signaux sources inconnus statistiquement indépendants l'un de l'autre. Le temps étant toujours un multiple entier de la période d'échantillonnage T,, on conviendra de noter x(t) le signal à temps discret correspondant au signal continu échantillonné à la date tT.. De cette façon, on pourra supposer dorénavant et sans restreindre la généralité de la présentation que T@ vaut 1.
xj(t) = s,(t) + B2 s2(t-#2) + v1(t), (1)
x2(t) = s2(t) + B1 s1(t-#1) + v2(t), (2)
où si (t) sont des signaux sources inconnus statistiquement indépendants l'un de l'autre. Le temps étant toujours un multiple entier de la période d'échantillonnage T,, on conviendra de noter x(t) le signal à temps discret correspondant au signal continu échantillonné à la date tT.. De cette façon, on pourra supposer dorénavant et sans restreindre la généralité de la présentation que T@ vaut 1.
La quantité #1 désigne le retard (ou l'avance si il est négatif) du signal si sur le capteur j, j X i,. Les quantités Bi sont introduites pour rendre compte éventuellement d'un amortissement de propagation, inférieur ou égal à I par définition. Si un des deux Bi disons par exemple B, est supérieur à 1, cela voudra simplement dire que les deux sources sont du même côté de la médiatrice à l'antenne. Dans ce cas, on pourra toujours poser B'2 = 1/B2' #'2 = -#2, et S 2(t) = B2 S2(t-T de sorte que:
x,(t) = s,(t) + s'2(t) + v,(t), (3) x1(t) = B2' s2'(t-#2') + B1 s1(t-#1) + v2(t), (4)
et que B'2 < 1.Enfin les termes vjt) sont des bruits Gaussiens inconnus éventuellement colorés et corrélés entre eux.
x,(t) = s,(t) + s'2(t) + v,(t), (3) x1(t) = B2' s2'(t-#2') + B1 s1(t-#1) + v2(t), (4)
et que B'2 < 1.Enfin les termes vjt) sont des bruits Gaussiens inconnus éventuellement colorés et corrélés entre eux.
Un modèle tout à fait équivalent aux précédents est le suivant:
x1(t) = s1(t) + B2 s2(t) + v1(t)
x2(t) = B, s1(t-#1) + s2(t-X2) + v2(t)
Le problème est donc le suivant: connaissant seulement un certain nombre de mesures {x1 (t), x2(t) i < t < T} , on veut estimer les amortissements B, et les retards x1, pour éventuellement en déduire les signaux sources si(t). On suppose que l'on ne connaît rien des propriétés spectrales des signaux sources, qui en outre peuvent fort bien être présents dans la même bande de fréquence.
x1(t) = s1(t) + B2 s2(t) + v1(t)
x2(t) = B, s1(t-#1) + s2(t-X2) + v2(t)
Le problème est donc le suivant: connaissant seulement un certain nombre de mesures {x1 (t), x2(t) i < t < T} , on veut estimer les amortissements B, et les retards x1, pour éventuellement en déduire les signaux sources si(t). On suppose que l'on ne connaît rien des propriétés spectrales des signaux sources, qui en outre peuvent fort bien être présents dans la même bande de fréquence.
Pour résoudre ce problème, I'invention propose un procédé d'estimation aveugle de retards différentiels entre deux signaux provenant de sources distinctes et reçus simultanément, principalement caractérisé en ce que l'on effectue une transformation de Fourier sur chacun de ces signaux, que l'on détermine les multispectres d'ordre 4 g1 à g6 pour une fréquence f des signaux ainsi transformés et les multispectres g1' à 6' pour la fréquence 2f1 que l'on détermine les deux premiers vecteurs:
a2 = B1 B2 g6 et: :
α2' = B1B2g6' que l'on résoud (110) le polynôme
α2' = B1B2g6' que l'on résoud (110) le polynôme
en retenant les deux racines dont le module est le plus proche de l'unité, et que l'on détermine les retards #1 et #2 des deux signaux reçus en résolvant les formules::
z1 = ej2#f#1, et z2 = ej2#f#2
Selon une autre caractéristique, on détermine si les deux deuxièmes vecteurs (g1, 92, 93) et (g'1, 9'2, 9'3) sont colinéaires et l'on résoud dans l'affirmative la formule
puis on calcule les coefficients d'amortissement Bi = #P et B2 =1/B1, et dans la négative on résoud les deux polynômes
(g2g1'-g2'g1)x+(-g3g1'-g3'g2)x+(g3g2'-g3'g2)=0, (g2 g'3 -g'2 g3)y2 +(- g1 g3-g1g3)y+(g1g2-g1g2)=0.
z1 = ej2#f#1, et z2 = ej2#f#2
Selon une autre caractéristique, on détermine si les deux deuxièmes vecteurs (g1, 92, 93) et (g'1, 9'2, 9'3) sont colinéaires et l'on résoud dans l'affirmative la formule
puis on calcule les coefficients d'amortissement Bi = #P et B2 =1/B1, et dans la négative on résoud les deux polynômes
(g2g1'-g2'g1)x+(-g3g1'-g3'g2)x+(g3g2'-g3'g2)=0, (g2 g'3 -g'2 g3)y2 +(- g1 g3-g1g3)y+(g1g2-g1g2)=0.
en ne retenant que les racines positives réelles et on calcule alors les coefficients d'amortissement Bi = W et B2 =
Selon une autre caractéristique, on détermine si le produit
B1 x B2 est différent de 1 et dans l'affirmative on résoud au sens des moindres carrés le système
pour déterminer finalement les retards entre les signaux à partir de la formule:#k=ej2#f#k.
Selon une autre caractéristique, on détermine si le produit
B1 x B2 est différent de 1 et dans l'affirmative on résoud au sens des moindres carrés le système
pour déterminer finalement les retards entre les signaux à partir de la formule:#k=ej2#f#k.
Selon une autre caractéristique, si le produit B1 x B2 est proche de l'unité on détermine si les deux vecteurs déterminés par
2a0 <SEP> 2a1
<tb> # <SEP> #z=# <SEP> #
<tb> <SEP> a1 <SEP> a2
<tb> sont colinéaires et dans l'affirmative on calcule leur facteur de proportionnalité z à partir de la formule ci-dessus pour déterminer finalement les retards entre les signaux à partir de la formule : Z = ej2#f#1 = ej2#f#2.
<tb> # <SEP> #z=# <SEP> #
<tb> <SEP> a1 <SEP> a2
<tb> sont colinéaires et dans l'affirmative on calcule leur facteur de proportionnalité z à partir de la formule ci-dessus pour déterminer finalement les retards entre les signaux à partir de la formule : Z = ej2#f#1 = ej2#f#2.
Selon une autre caractéristique, si ces deux premiers vecteurs ne sont pas colinéaires, on détermine si les deux quatrièmes vecteurs déterminés par
α0 <SEP> alpha;2
<tb> sont colinéaires et dans l'affirmative on calcule leur facteur de proportionnalité z à partir de la formule ci-dessus pour déterminer finalement les retards entre les signaux par la formule #z=ej2#f#1=-ej2#f#2.
<tb> sont colinéaires et dans l'affirmative on calcule leur facteur de proportionnalité z à partir de la formule ci-dessus pour déterminer finalement les retards entre les signaux par la formule #z=ej2#f#1=-ej2#f#2.
Selon une autre caractéristique, les transformés de Fourier et la détermination des multispectres ont effectuées de manière récursive.
Selon une autre caractéristique, on reçoit les signaux sur une antenne linéaire acoustique comprenant un ensemble d'hydrophones, on forme une voie "end-fire" et on applique le procédé selon l'une quelconque des revendications 1 à 6 entre le signal de la voie "end-fire" et chacun des signaux reçus par les hydrophones de l'antenne.
D'autres particularités et avantages de l'invention apparaîtront clairement dans la description suivante, présentée à titre d'exemple non limitatif en regard des figures annexées qui représentent:
- la figure 1, un schéma synoptique général d'un dispositif mettant en oeuvre l'invention;
- la figure 2, un mode de réalisation de l'unité 103 de la figure 1;
- la figure 3, un mode de réalisation par calcul récursif; et
- la figure 4, un exemple d'application de l'invention en acoustique sous-marine.
- la figure 1, un schéma synoptique général d'un dispositif mettant en oeuvre l'invention;
- la figure 2, un mode de réalisation de l'unité 103 de la figure 1;
- la figure 3, un mode de réalisation par calcul récursif; et
- la figure 4, un exemple d'application de l'invention en acoustique sous-marine.
Le dispositif représenté en figure 1, comporte deux unités de calcul de Transformée de Fourier 101 et 102, une unité de calcul des multispectres 103 (spectres cumulants), que nous appelerons unité de calcul préliminaire, deux unités de calcul des amortissements 104 (A1) et 105 (A2), et cinq unités de calcul des retards 106 (R1) à 110 (R5), deux unités de test 111 (T1) et 112 (T2), et un ensemble de commutateurs 113 à 117.
Le contenu et le fonctionnement de ces organes sont décrits ciaprès.
Les signaux reçus xl (t) et x2 (t) sont donc traités dans les unités de calcul de transformée de Fourier rapide FFT 101 et 102 pour obtenir des signaux exprimés en fréquences X1 (f) et X2 (f).
Ils sont ensuite appliqués à l'unité de calcul des multispectres d'ordre 4 g1 à g6 suivants (on trouvera dans D.R. BRILLINGER "Time Series,
Data Analysis and Theory" Holden-Day, 1981 [1] la définition de ces multispectres):
g1 = C{X1(f),X1(f)*,X1(f),X1(f)*} (5)
g2 = C{X1(f),X1(f)*,X2(f),X2(f)*} (6)
g3 = C{X2(f),X2(f)*,X2(f),X2(f)*} (7)
g4 = C{X1(f),X1(f)*,X1(f),X2(f)*} (8)
g, = C{X1(f),X2(f)*,X2(f),X2(f)*} (9)
g6 = C{X1(f),X2(f)*,X1(f),X2(f)*} (10)
où (*) désigne la conjugaison complexe, Xi(f) la transformée de
Fourier de xi(t) obtenue dans 101 ou 10?, et f la fréquence dans
Data Analysis and Theory" Holden-Day, 1981 [1] la définition de ces multispectres):
g1 = C{X1(f),X1(f)*,X1(f),X1(f)*} (5)
g2 = C{X1(f),X1(f)*,X2(f),X2(f)*} (6)
g3 = C{X2(f),X2(f)*,X2(f),X2(f)*} (7)
g4 = C{X1(f),X1(f)*,X1(f),X2(f)*} (8)
g, = C{X1(f),X2(f)*,X2(f),X2(f)*} (9)
g6 = C{X1(f),X2(f)*,X1(f),X2(f)*} (10)
où (*) désigne la conjugaison complexe, Xi(f) la transformée de
Fourier de xi(t) obtenue dans 101 ou 10?, et f la fréquence dans
<SEP> 1 <SEP> 1
<tb> 2 > 2 <SEP>
<tb>
On définit également les mêmes quantités pour la fréquence double à partir de X,{2f), que l'on note g:, g2' et g3' respectivement. Le signe C désigne le cumulant conjoint des variables figurant en argument. La figure 2 illustre l'implantation du calcul des g1 pour une fréquence arbitraire f. C'est le même dispositif qui est utilisé pour la fréquence 2f.
<tb> 2 > 2 <SEP>
<tb>
On définit également les mêmes quantités pour la fréquence double à partir de X,{2f), que l'on note g:, g2' et g3' respectivement. Le signe C désigne le cumulant conjoint des variables figurant en argument. La figure 2 illustre l'implantation du calcul des g1 pour une fréquence arbitraire f. C'est le même dispositif qui est utilisé pour la fréquence 2f.
On décrira plus loin dans ce texte une autre implantation de nature récursive.
Dans la figure 2, le signal X1 est appliqué à un multiplicateur 201 et à un organe de calcul du module carré 202. Le signal X2 est appliqué au multiplicateur 201 par l'intermédiaire d'un organe de calcul du complexe conjugué 203 et à un autre organe de calcul 204 identique à l'organe 202.
La sortie de l'organe 202 est appliquée à un multiplicateur 205, à un organe de calcul du carré 206 et à un organe de moyennage statistique 207.
La sortie du multiplicateur 201 est appliquée au multiplicateur 205, à un organe 208 identique à l'organe 202, à un organe 209 identique à l'organe 206, à un multiplicateur 210, et à un organe 211 identique à l'organe 207.
La sortie de l'organe 204 est appliquée au multiplicateur 240, à un organe 212 identique à l'organe 206 et à un organe 213 identique à l'organe 207.
Les sorties des multiplicateurs 205 et 210 et celles des organes 206, 208, 209, 212 sont respectivement appliquées à des organes 214 à 219 identiques à l'organe 207.
La sortie de l'organe 207 est appliquée à deux multiplicateurs 220 et 221 ainsi qu'à un organe 222 identique à l'organe 206.
La sortie de l'organe 213 est reliée au multiplicateur 221, à un autre multiplicateur 223, et à un organe 224 identique à l'organe 202.
La sortie de l'organe 211 est reliée aux multiplicateurs 220 et 223, à un organe 225 identique à l'organe 206 et à un organe 226 identique à l'organe 202.
Les sorties de l'organe 225 et du multiplicateur 221 sont reliées respectivement à des multiplicateurs scalaires par -1, 227 et 228.
Les sorties des multiplicateurs 220 et 223 et les organes 222, 224, 225, 226 sont reliées respectivement à des multiplicateurs scalaires par -2, 229 à 233.
Les sorties des organes 214 à 219 sont reliées respectivement à des additionneurs 234 à 239 qui sont aussi reliés respectivement aux sorties des organes 228 à 233.
La sortie de l'additionneur 236 est reliée à un additionneur 240 qui est aussi relié à la sortie de l'organe 227. Ces additionneurs 234 à 240, à l'exception du 236, délivrent les spectres cumulants g1 à 96.
L'estimation des facteurs d'amortissement se fait de deux manières différentes, selon que les vecteurs des multispectres de fréquence f et 2f sont colinéaires ou non.
Si le vecteur (g1, g2, g3) est colinéaire à (g1',g2',g3'), alors on calcule dans l'unité de calcul Al
et on déduit B1 = #p et B2 = 1/B1. Le test de colinéarité peut se faire par exemple par l'intermédiaire du cosinus de l'angle, rapport du produit scalaire des deux vecteurs sur le produit de leur norme, dans l'unité de test 111. Selon les résultats le vecteur est dirigé vers l'unité A1 ou l'unité
A2 pour un commutateur 113..
A2 pour un commutateur 113..
Sinon, on résoud dans l'unité de calcul A2 alimentée par le commutateur A1 les deux polynômes du second degré suivants:
(g2 g'1 -g'2 gl )x2 +(-g3g1' -g'3 g2 )x + (g3 g'2 -g'3 g2)=0, (12) (g2g3' -g'2 g3)y2 + (-g1g3' 3 g3)y+(g1g2' -g'1 g2) = o. (13)
On ne retiendra que les racines positives réelles, et on en déduira
B1 = #x et B2 = #y.
(g2 g'1 -g'2 gl )x2 +(-g3g1' -g'3 g2 )x + (g3 g'2 -g'3 g2)=0, (12) (g2g3' -g'2 g3)y2 + (-g1g3' 3 g3)y+(g1g2' -g'1 g2) = o. (13)
On ne retiendra que les racines positives réelles, et on en déduira
B1 = #x et B2 = #y.
Les amortissements B, et B2 sont ensuite utilisés pour calculer les retards #i. Toutefois les valeurs de ces paramètres sont disponibles sur les sorties des unités A1 et A2, ce qui permet de vérifier directement les valeurs obtenues, par exemple pour effectuer des tests de vraisemblance.
Un commutateur 114 fonctionnant en phase avec le commutateur 113 permet de sélectionner la sortie de A1 ou de A2 selon le calcul effectué en fonction de la colinéarité des vecteurs.
Pour éliminer l'effet des facteurs d'amortissement B1 et B2 dans les calculs, on effectue dans l'unité de test 112 le calcul du produit B1 x B2.
Si ce produit B,x B2 est différent de 1 alors on envoie les vecteurs dans l'unité de calcul R1 par les commutateurs 115 et les retards sont calculés directement par l'unité R1 ou R6. On résoud au sens des moindres carrés le système linéaire surdéterminé suivant, toujours de rang plein puisque B1 x B2 $ 1:
Ceci permet d'obtenir les multispectres des sources à la fréquence f, notés yi et γ2 ici.Puis on résoud, toujours dans R1, le système suivant dans le corps des complexes:
<tb> γ1 <SEP> <SEP> B1 <SEP> γ2 <SEP> B2 <SEP> #1 <SEP> <SEP> g4
<tb> # <SEP> ## <SEP> #=# <SEP> # <SEP> (15)
<tb> <SEP> γ1 <SEP> B1 <SEP> γ2 <SEP> B2 <SEP> #2 <SEP> g5
<tb>
Les retards nc sont alors obtenus en résolvant la relation: = = eJ (16)
Toutes ces opérations sont effectuées dans R1, qui délivre directement #1 et #2.
<tb> # <SEP> ## <SEP> #=# <SEP> # <SEP> (15)
<tb> <SEP> γ1 <SEP> B1 <SEP> γ2 <SEP> B2 <SEP> #2 <SEP> g5
<tb>
Les retards nc sont alors obtenus en résolvant la relation: = = eJ (16)
Toutes ces opérations sont effectuées dans R1, qui délivre directement #1 et #2.
Sinon, lorsque le produit B, B2 est proche de l'unité, on envoie les vecteurs dans l'unité de calcul R2 par le commutateur 115 et l'on calcule d'abord les quantités suivantes pour la fréquence f
α2 = B1B2g6 (19) et les mêmes quantités pour la fréquence 2f.
α2 = B1B2g6 (19) et les mêmes quantités pour la fréquence 2f.
a'2 = B1B2g'6 (22)
Les résultats obtenus dans l'unité R2 sont ensuite utilisés pour estimer les retards dans les unités R3, R4 ou R5.
Les résultats obtenus dans l'unité R2 sont ensuite utilisés pour estimer les retards dans les unités R3, R4 ou R5.
Si les vecteurs définis cidessous sont colinéaires, on calcule alors leur facteur de proportionnalité z:
<tb> #2a0# <SEP> #2a1#
<tb>
Puis on calcule les retards par:
Z = ej2#f#1 = eJ2 (24)
L'estimation de la colinéarité et les calculs du facteur de proportionnalité sont effectués dans l'unité R3 qui délivre #1 et X2 par l'intermédiaire d'un commutateur 116 actionné par R3 selon le résultat du test de colinéarité.
<tb>
Puis on calcule les retards par:
Z = ej2#f#1 = eJ2 (24)
L'estimation de la colinéarité et les calculs du facteur de proportionnalité sont effectués dans l'unité R3 qui délivre #1 et X2 par l'intermédiaire d'un commutateur 116 actionné par R3 selon le résultat du test de colinéarité.
Dans le cas où les vecteurs ci-dessus ne sont pas colinéaires,on cherche si les vecteurs ci-dessous sont colinéaires, et on calcule alors leur coefficient de proportionnalité z:
<SEP> α0 <SEP> <SEP> α2
<tb> #α0'#z=#α1'# <SEP> (25)
<tb> α1' <SEP> α2'
<tb>
Puis on calcule les retards par:
= & nft1 = -ej2#f#2 (26)
Ces calculs sont effectués dans l'unité R4 qui délivre #1 et #2 par l'intermédiaire d'un commutateur 117 actionné par R4 selon le résultat du test de colinéarité.
<tb> #α0'#z=#α1'# <SEP> (25)
<tb> α1' <SEP> α2'
<tb>
Puis on calcule les retards par:
= & nft1 = -ej2#f#2 (26)
Ces calculs sont effectués dans l'unité R4 qui délivre #1 et #2 par l'intermédiaire d'un commutateur 117 actionné par R4 selon le résultat du test de colinéarité.
Enfin, si aucun des deux cas précédents n'a été rencontré, on résoud le polynôme de degré 4 suivant: (a,2a'1 -a,2aO)z4 +2a1(α2α0'-α0α1')z +(2al2a', -aO2a'2 -a22a'0)z2 (27) +2a1(α0α2'-α2α1')z + (a22a', a12a'2) = O (28)
Parmi les 4 racines de ce polynôme, on retiend les deux dont le module est le plus proche de l'unité, et on obtient les valeurs des retards par:
z1 = ej2#f#1, et z2 = ej2#f#2 (29)
Ces calculs sont effectués dans l'unité R5 qui peut traiter la plupart des cas de figure, les unités R3 et R4 ne traitant que deux cas particuliers. Dans un mode de réalisation simplifié, on pourra donc se passer des unités R3 et R4. En outre, si les amortissements peuvent être négligés, on pourra aussi se passer des unités A1 et A2. Dans un tel dispositif simplifié, les unités à conserver impérativement seront donc les unités R2 et
R5, en prenant B1 = B2 = 1.
Parmi les 4 racines de ce polynôme, on retiend les deux dont le module est le plus proche de l'unité, et on obtient les valeurs des retards par:
z1 = ej2#f#1, et z2 = ej2#f#2 (29)
Ces calculs sont effectués dans l'unité R5 qui peut traiter la plupart des cas de figure, les unités R3 et R4 ne traitant que deux cas particuliers. Dans un mode de réalisation simplifié, on pourra donc se passer des unités R3 et R4. En outre, si les amortissements peuvent être négligés, on pourra aussi se passer des unités A1 et A2. Dans un tel dispositif simplifié, les unités à conserver impérativement seront donc les unités R2 et
R5, en prenant B1 = B2 = 1.
Par ailleurs, la transformée de Fourier d'une part, et les multispectres d'autre part, peuvent être mis à jour récursivement à chaque fois qu'un nouvel échantillon [x,(t)x2(t)] est reçu. Ainsi si on le désire,
L'ensemble de l'unité de calcul préliminaire peut fonctionner de façon récursive.
L'ensemble de l'unité de calcul préliminaire peut fonctionner de façon récursive.
Ce calcul récursif est surtout intéressant lorsque seules certaines bandes de fréquence sont utiles, ou bien lorsque le signal ne peut être mémorisé sur la durée complète d'une transformée de Fourier. En notant
X(f) la TF d'un signal x(t) dont les échantillons sont successivement observés, si N est la longueur de la TF, on a à une date t quelconque:
si on convient de dire que la fréquence varie dans l'ensemble
X(f) la TF d'un signal x(t) dont les échantillons sont successivement observés, si N est la longueur de la TF, on a à une date t quelconque:
si on convient de dire que la fréquence varie dans l'ensemble
Le calcul récursif proposé est:
X(f;t+1)={ss[X(f;t)-x(t-N+1)]+x(t+1)}e2#jf (31)
où ss est un nombre positif très proche de 1 dont le but est d'éviter les explosions numériques.
X(f;t+1)={ss[X(f;t)-x(t-N+1)]+x(t+1)}e2#jf (31)
où ss est un nombre positif très proche de 1 dont le but est d'éviter les explosions numériques.
Par exemple, sur une machine en simple précision, on pourra prendre ss= l - 104.
Si l'ensemble des N échantillons {x(t-N+1),..,x(t)} est stocké dans un registre à décalage comme représenté en figure 3, I'échantillon x(t-N+1) sera directement disponible pour le calcul de X(f;t + 1) avant d'être définitivement oublié.
Sur la figure 3, on a représenté un dispositif permettant de faire ce calcul avec un registre à décalage 301 recevant x (t + 1). La sortie est multipliée par le facteur ss dans un multiplicateur 302. Un additionneur 303 reçoit la sortie de 302, x(t+1) et la sortie X(f) retardée de T dans 304 et multipliée par ss dans 305. Une unité de calcul par e2xjf 306 est connectée en sortie de 303 et délivre x(f) sur la sortie et à 304.
Cependant, une autre variante est éventuellement possible, et présente un intérêt si tous les échantillons spectraux X(f;t) sont calculés.
Dans ce cas, au lieu de mémoriser ces N échantillons, on peut recalculer x(t-N+1) à partir de X(f;t) comme suit:
Ceci procure un gain en place mémoire (il n'est plus nécessaire d'avoir recours à un registre à décalage), au prix d'une augmentation de coût calcul. La stabilité numérique est encore assurée par la présence du facteur de fuite ss.
Pour des observations centrées, le calcul récursif des multispectres s'effectue alors en calculant d'une part les spectres à la fréquence f à l'aide des récurrences: 11(t+1)=(1-#) 11(t)+#|X1(f;t)| (33)
12(t+1)=(1-#) 12(t)+#(X1(f;t)X2(f;t)*) (34)
22(t+1)=(1-#) 22(t)+#|X2(f;t)| (35)
Le coefficient E est nombre réel de l'intervalle ]0,1[ déterminant la puissance du moyennage (et donc la longueur de la mémoire); il est parfois appelé "facteur d'oubli".D'autre part on calcule les multispectres à la fréquence f, gk, avec les formules: g1(t+1)=(1-#)g1(t)+#(|X1(f;t)|4-2 11(t)) (36) g2(t+1)=(1-#)g2(t)+#(|X1(f;t)||X2(f;t)|- 11(t) 22(t)-| 12(t)|) (37) g3(t+1)=(1-#)g3(t)+#(|X2(f;t)|4-2 22(t)) (38) g4(t + 1) = (1 - #)g4(t)+#(|X1(f;t)|X1(f;t)X2(f;t)*-2 11(t) 12(t)) (39) g5(t + 1) = (1 - #)g5(t) + rkxî(j;t)x2(j;t) *jX2(f;t)J2 -2 12(t) 22(t)) (40) g6(t+1)=(1-#)g6(t)+#(X1(f;t)X2(f;t)*-2 12(t)) (41)
Ici encore, le coefficient n est un nombre réel de l'intervalle ]0,1[ déterminant la puissance du moyennage.Il peut éventuellement être égal à
E. Les multispectres à la fréquence double, g'k, sont définis de la même manière.
12(t+1)=(1-#) 12(t)+#(X1(f;t)X2(f;t)*) (34)
22(t+1)=(1-#) 22(t)+#|X2(f;t)| (35)
Le coefficient E est nombre réel de l'intervalle ]0,1[ déterminant la puissance du moyennage (et donc la longueur de la mémoire); il est parfois appelé "facteur d'oubli".D'autre part on calcule les multispectres à la fréquence f, gk, avec les formules: g1(t+1)=(1-#)g1(t)+#(|X1(f;t)|4-2 11(t)) (36) g2(t+1)=(1-#)g2(t)+#(|X1(f;t)||X2(f;t)|- 11(t) 22(t)-| 12(t)|) (37) g3(t+1)=(1-#)g3(t)+#(|X2(f;t)|4-2 22(t)) (38) g4(t + 1) = (1 - #)g4(t)+#(|X1(f;t)|X1(f;t)X2(f;t)*-2 11(t) 12(t)) (39) g5(t + 1) = (1 - #)g5(t) + rkxî(j;t)x2(j;t) *jX2(f;t)J2 -2 12(t) 22(t)) (40) g6(t+1)=(1-#)g6(t)+#(X1(f;t)X2(f;t)*-2 12(t)) (41)
Ici encore, le coefficient n est un nombre réel de l'intervalle ]0,1[ déterminant la puissance du moyennage.Il peut éventuellement être égal à
E. Les multispectres à la fréquence double, g'k, sont définis de la même manière.
Ce procédé permet donc de séparer les signaux de deux sources de bruit distinctes si celles-ci sont statistiquement indépendantes. C'est plus particulièrement le cas en acoustique sous-marine pour le rayonnement acoustique de deux bateaux, mais il existe de nombreuses autres situations où cette condition est remplie, par exemple pour deux avions.
En outre ce procédé est insensible au bruit Gaussien, qui est éliminé, ce qui est tout à fait précieux pour repérer des sources dont le bruit propre est noyé dans le bruit ambiant, comme c'est aussi souvent le cas en acoustique sous-marine.
Une application particulière représentée en figure 4 permet d'éliminer le bruit du bateau tracteur 403 dans le cas d'une antenne acoustique linéaire remorquée 401. Pour cela on forme tout d'abord une voie 404 dite "end-fire" dont le lobe principal est dirigé selon l'axe de l'antenne vers le bateau tracteur, en sommant les signaux des hydrophones 402 retardés chacun de d/c, où d est la distance entre ces hydrophone, et c la vitesse du son dans l'eau. On applique ensuite le procédé de l'invention entre le signal de cette voie et chaque signal provenant des hydrophones.
Ceci permet de séparer les signaux provenant du bateau tracteur de ceux provenant d'une autre source de bruit, tel qu'un autre bateau que l'on cherche à repérer. En outre on élimine ainsi le bruit sous-marin ambiant.
On peut aussi utiliser le procédé de l'invention pour réduire le bruit dans une communication acoustique entre un bateau de surface et un véhicule sous-marin, ou pour discriminer les émissions de deux sonars actifs.
Enfin le procédé selon l'invention peut aussi s'utiliser dans les techniques radar, d'une manière semblable à celle utilisée dans les techniques sonar.
Claims (7)
1. Procédé d'estimation aveugle de retards différentiels entre deux signaux provenant de sources distinctes et reçus simultanément, caractérisé en ce que l'on effectue une transformation de Fourier (101, 102) sur chacun de ces signaux, que l'on détermine (103) les multispectres d'ordre 4 g1 à g6 pour une fréquence f des signaux ainsi transformés et les multispectres g1' à g6' pour la fréquence 2f1, que l'on détermine (107) les deux premiers vecteurs:
α2=B1B2g6 et::
α2'=B1B3g6' que l'on résoud (110) le polynôme
(g2g3'-g2'g3)y+(-g1g3'-g1'g3)+(g1g2'-g1'g2)=0.
(g2g1'-g2'g1)x+(-g3g1'-g3'g2)x+(g3g2'-g3'g2)=0,
puis que l'on calcule les coefficients d'amortissement B1 = #p et B2 =1/B1, et que dans la négative on résoud (105) les deux polynômes
2.Procédé selon la revendication 1, caractérisé en ce que l'on détermine (111) si les deux deuxièmes vecteurs (g1, 92, 93) et g1', 9'2, 9'3) sont colinéaires et que l'on résoud (104) dans l'affirmative la formule
en retenant les deux racines dont le module est le plus proche de l'unité, et que l'on détermine (110) les retards 71 et #2 des deux signaux reçus en résolvant les formules: z1=ej2#f#1,et z2=ej2#f#2
en ne retenant que les racines positives réelles et que l'on calcule alors les coefficients d'amortissement Bi = W et B2 =
3. Procédé selon la revendication 2, caractérisé en ce que l'on détermine si le produit B1 x B2 est différent de 1 et que dans l'affirmative on résoud (112) au sens des moindres carrés le système
4.Procédé selon la revendication 3, caractérisé en ce que si le produit B1 x B2 est proche de l'unité on détermine (108) si les deux vecteurs déterminés par
pour déterminer finalement les retards entre les signaux à partir de la formule:#k=
<tb> sont colinéaires et que dans l'affirmative on calcule (108) leur facteur de proportionnalité z à partir de la formule cidessus pour déterminer finalement (108) les retards entre les signaux à partir de la formule: Z = ej2+rl = e'l.
<tb> <SEP> a1 <SEP> a2
<tb> # <SEP> #z=# <SEP> #
<tb> 2a0 <SEP> 2a1
5. Procédé selon la revendication 4, caractérisé en ce que si ces deux premiers vecteurs ne sont pas colinéaires, on détermine (109) si les deux quatrièmes vecteurs déterminés par
<tb> sont colinéaires et que dans l'affirmative on calcule (109) leur facteur de proportionnalité z à partir de la formule ci-dessus pour déterminer finalement (109) les retards entre les signaux par la formule #z=ej2#f#1=-ej2#f#2.
<tb> <SEP> α1' <SEP> <SEP> α2'
<tb> <SEP> α0 <SEP> <SEP> α2
6. Procédé selon l'une quelconque des revendications 1 à 5, caractérisé en ce que les transformés de Fourier (101, 102) et la détermination des multispectres (103) sont effectuées de manière récursive.
7. Procédé selon l'une quelconque des revendications 1 à 6, caractérisé en ce que l'on reçoit les signaux sur une antenne linéaire (101) acoustique comprenant un ensemble d'hydrophones, que l'on forme une voie "end-fire" (104) et que l'on applique le procédé selon l'une quelconque des revendications 1 à 6 entre le signal de la voie "end-fire" et chacun des signaux reçus par les hydrophones (102) de l'antenne.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| FR9410325A FR2724028B1 (fr) | 1994-08-26 | 1994-08-26 | Procede d'estimation aveugle de retards differentiels entre deux signaux |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| FR9410325A FR2724028B1 (fr) | 1994-08-26 | 1994-08-26 | Procede d'estimation aveugle de retards differentiels entre deux signaux |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| FR2724028A1 true FR2724028A1 (fr) | 1996-03-01 |
| FR2724028B1 FR2724028B1 (fr) | 1996-09-20 |
Family
ID=9466510
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| FR9410325A Expired - Fee Related FR2724028B1 (fr) | 1994-08-26 | 1994-08-26 | Procede d'estimation aveugle de retards differentiels entre deux signaux |
Country Status (1)
| Country | Link |
|---|---|
| FR (1) | FR2724028B1 (fr) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP0152523A2 (fr) * | 1984-02-22 | 1985-08-28 | Messerschmitt-Bölkow-Blohm Gesellschaft mit beschränkter Haftung | Procédé pour localiser des émetteurs de signaux avec suppression des interférences |
| JPH02115780A (ja) * | 1988-10-25 | 1990-04-27 | Nec Corp | 送波レベル測定装置 |
| JPH04335175A (ja) * | 1991-05-10 | 1992-11-24 | Nec Corp | 目標検出装置 |
| US5323459A (en) * | 1992-11-10 | 1994-06-21 | Nec Corporation | Multi-channel echo canceler |
-
1994
- 1994-08-26 FR FR9410325A patent/FR2724028B1/fr not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP0152523A2 (fr) * | 1984-02-22 | 1985-08-28 | Messerschmitt-Bölkow-Blohm Gesellschaft mit beschränkter Haftung | Procédé pour localiser des émetteurs de signaux avec suppression des interférences |
| JPH02115780A (ja) * | 1988-10-25 | 1990-04-27 | Nec Corp | 送波レベル測定装置 |
| JPH04335175A (ja) * | 1991-05-10 | 1992-11-24 | Nec Corp | 目標検出装置 |
| US5323459A (en) * | 1992-11-10 | 1994-06-21 | Nec Corporation | Multi-channel echo canceler |
Non-Patent Citations (5)
| Title |
|---|
| HSING HSING CHIANG ET CHRYSOSTOMOS L. NIKIAS: "A new method for adaptative time delay estimation for non gaussian signals", IEEE TRANSACTIONS ON ACOUSTICS , SPEECH AND SIGNAL PROCESSING, vol. 38, no. 2, February 1990 (1990-02-01), NEW YORK , NY , USA ., pages 209 - 19, XP000113160 * |
| HSING HSING CHIANG ET CHRYSOSTOMOS L. NIKIAS: "Cumulant based adaptative time delay estimation", CONFERENCE RECORD OF THE TWENTY SECOND ASILOMAR CONFERENCE ON SIGNAL , SYSTEMS AND COMPUTERS, vol. 1, 31 October 1988 (1988-10-31), PACIFIC GROVE , CA , USA, pages 15 - 9, XP000130213 * |
| JITENDRA K. TUGNAIT: "On time delay estimation with unknown spatially correlated gaussian noise using fourth order cumulants and cross cumulants", IEEE TRANSACTIONS ON SIGNAL PROCESSING, vol. 39, no. 6, June 1991 (1991-06-01), NEW YORK , NY , USA, pages 1258 - 67 * |
| PATENT ABSTRACTS OF JAPAN vol. 14, no. 337 (P - 1079) 20 July 1990 (1990-07-20) * |
| PATENT ABSTRACTS OF JAPAN vol. 17, no. 182 (P - 1518) 8 April 1993 (1993-04-08) * |
Also Published As
| Publication number | Publication date |
|---|---|
| FR2724028B1 (fr) | 1996-09-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| EP0605281B1 (fr) | Procéde de débruitage vectoriel de la parole et dispositif de mise en oeuvre | |
| EP0665665B1 (fr) | Procédé et dispositif permettant à un modem de se synchroniser sur un transmetteur de données numériques par voie hertzienne en présence de brouilleurs | |
| EP2538409B1 (fr) | Procédé de débruitage pour équipement audio multi-microphones, notamment pour un système de téléphonie "mains libres" | |
| EP2100161B1 (fr) | Procede de traitement radar passif multivoies d'un signal d'opportunite en fm | |
| EP0867731B1 (fr) | Procédé et dispositif de levée d'ambiguité en distance appliquée notamment à un radar à saut de fréquence | |
| CA2400763A1 (fr) | Procede de localisation de sources radioelectriques au moyen d'un radiogoniometre haute resolution deux voies | |
| EP1615049B1 (fr) | Traitement cohérent rapide pour codes à spectre de raies périodiques | |
| EP0749083A1 (fr) | Procédé et dispositif de détermination du spectre de fréquence d'un signal | |
| CA1206584A (fr) | Dispositif de filtrage adaptatif de signaux recus par un sonar actif pour la rejection de la reverberation | |
| EP0068909B1 (fr) | Procédé et dispositif de réduction de la puissance des signaux de brouillage reçus par les lobes secondaires d'une antenne radar | |
| FR2483086A1 (fr) | Procede de traitement de signal pour radar a visee laterale et a synthese d'ouverture et circuit de mise en oeuvre | |
| FR2584213A1 (fr) | Dispositif de calcul d'une transformee de fourier discrete, glissante, et son application a un systeme radar. | |
| FR2652654A1 (fr) | Echographe ultrasonore utilisant un dispositif numerique de formation de voies en reception. | |
| FR2773283A1 (fr) | Dispositif et procede de traitement numerique de signaux | |
| CA1226654A (fr) | Dispositif de discrimination d'echos radar | |
| EP0334711B1 (fr) | Dispositif d'élimination du fouillis mobile dans un radar | |
| EP0464115A1 (fr) | Procede et dispositif d'analyse spectrale en temps reel de signaux instationnaires complexes | |
| FR2724028A1 (fr) | Procede d'estimation aveugle de retards differentiels entre deux signaux | |
| FR2761550A1 (fr) | Filtre numerique pour retards fractionnaires | |
| EP1155497B1 (fr) | Procede et systeme de traitement de signaux d'antenne | |
| EP0333566B1 (fr) | Procédé et dispositifs de translation le long de l'axe des fréquences du module de la fonction de transfert d'un filtre | |
| CA2031719C (fr) | Procede et dispositif pour imposer un diagramme de rayonnement au repos a un reseau d'antennes de reception a formation adaptative de faisceau par le calcul | |
| FR2599855A1 (fr) | Dispositif de levee d'ambiguite de distance et de vitesse plus particulierement utilise associe a un recepteur radar du type doppler, et radar le comportant | |
| EP0029760B1 (fr) | Dispositif de filtrage et son application aux systèmes de visualisation radar | |
| FR2678738A1 (fr) | Procede de detection autoregressive d'un signal sinusouidal complexe dans du bruit et d'estimation de sa frequence pour un radar a impulsions. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| ST | Notification of lapse |