Brevet FR2926148A1 — publié le 20090710

PROCEDE DE SEPARATION DE SIGNAUX MELANGES EN UNE PLURALITE DE SIGNAUX COMPOSANTS La présente invention concerne le domaine technique général du traitement de signaux, et plus précisément de la séparation de signaux sources à partir de mélanges de signaux sources en utilisant la séparation aveugle de sources. La présente invention peut trouver de nombreuses applications. Elle peut s'appliquer notamment au traitement d'électroencéphalogrammes, au traitement de signaux dans le domaine des radiocommunications, ou encore au traitement de spectres de masse. PRESENTATION GENERALE DE L'ART ANTERIEUR Le problème de séparation de sources est un problème relativement récent en traitement du signal, qui consiste à séparer des sources, statistiquement indépendantes, à partir de plusieurs mélanges enregistrés soit par un réseau de capteurs à un instant donné (dispersion spatiale), soit par un capteur unique à différents instant (dispersion temporelle). Il arrive fréquemment qu'un phénomène physique observé soit la superposition de plusieurs phénomènes élémentaires indépendants et que cette superposition se traduise par une combinaison (ou mélange) linéaire instantanée. Les données de mesure issues de l'observation de ce phénomène sont alors un mélange linéaire instantané des contributions de chaque phénomène élémentaire indépendant. En plaçant une pluralité de capteurs, il est alors possible de mesurer plusieurs mélanges linéaires instantanés différents de ces contributions. Par exemple, en électro-encéphalographie, diverses électrodes sont placées sur le cuir chevelu d'un patient. Ces électrodes permettent de mesurer l'activité cérébrale du cerveau du patient. Chaque électrode mesure un mélange de signaux électriques émis par différentes sources du cerveau, les différentes électrodes mesurant différents mélanges. Des mélanges de sources sont dits différents lorsque les sources y apparaissent dans des proportions ou avec des pondérations différentes. Dans ce cas, aucune autre information n'est connue : - sur les sources et - sur la façon dont les signaux issus de ces sources se mélangent. On utilise alors la séparation de source pour identifier les pondérations du mélange dans un premier temps, et pour séparer dans les mesures la contribution de chaque source dans un second temps. Ceci permet une analyse séparée des différentes contributions sans être gêné par le mélange. Il existe actuellement deux types de méthodes pour la séparation de source. Le premier type de méthodes utilise une étape préalable dite de � blanchiment � permettant de faciliter les étapes de calcul ultérieur pour identifier les pondérations du mélange. Cependant, cette étape de blanchiment induit une erreur qui ne peut être corrigée lors de la séparation de la contribution de chaque source. Ce type de méthodes est donc peu précis. Le deuxième type de méthodes ne nécessite pas d'étape de blanchiment ; toutefois, ce deuxième type de méthode pose des problèmes en termes d'efficacité et de coût de calcul. Un but de la présente invention est de proposer un procédé de séparation de signaux mélangés en une pluralité de signaux composants permettant de pallier au moins l'un des inconvénients des types de méthodes précités. PRESENTATION DE L'INVENTION A cet effet on prévoit un procédé de séparation de signaux mélangés (xi(t), ..., xN(t)) enregistrés par au moins deux capteurs en une pluralité de signaux composants (si(t), ..., sM(t)), les signaux mélangés (xi(t), ..., xN(t)) correspondant à N combinaisons linéaires des M signaux composants (si(t), sM(t)) issus de sources indépendantes, chaque signal mélangé étant apte à s'écrire comme le produit d'une matrice de mélange (A) multiplié par un vecteur des signaux composants issus des sources (si(t), .., sM(t)), le procédé comprenant en outre une étape de diagonalisation conjointe dans laquelle la matrice de mélange A est estimée à partir d'un produit de matrices de rotation correspondant à une ou des matrices de rotation de Givens (G(8)) et à une ou des matrices de � rotation hyperbolique � (H(cp)). On entend, dans le cadre de la présente invention, par � matrice de rotation hyperbolique �, une matrice de rotation comportant des termes qui sont des fonctions hyperboliques (par exemple cosh, sinh) et dont le déterminant est égal à 1. Ainsi, l'invention permet d'identifier la matrice de mélange par une seule étape d'optimisation rapide et fiable. Le fait qu'il n'y ait qu'une seule étape fournit un résultat plus précis que les méthodes en deux étapes (blanchiment puis diagonalisation conjointe) telles que [1], [1b], [2]. En effet, les erreurs inévitables de la première étape ne peuvent être corrigées par la seconde étape. De plus, la méthode proposée utilise un algorithme algébrique (il s'agit de produits de rotations � trigonométriques � et � hyperboliques � de taille 2x2) particulièrement robuste et peu couteux en terme de puissance de calcul. Les méthodes comprenant une étape de blanchiment construisent une matrice de mélange orthogonale par produit de matrices orthogonales élémentaires (rotations de Givens). La présente invention permet quant à elle de construire une matrice de mélange de déterminant égal à 1 par multiplication de matrices de 3 déterminant égal à 1. Le groupe des matrices orthogonales (de déterminant positif +1) est inclus dans le groupe des matrices de déterminant égal à 1. Les étapes de calcul de la matrice de mélange sont donc facilitées avec la méthode de la présente invention qui permet de limiter les contraintes initiales sur les mélanges. L'étape initiale de � blanchiment � n'est plus indispensable mais la robustesse et l'efficacité des méthodes � algébriques � sont conservées. Par conséquent, la méthode de la présente invention permet le calcul d'une matrice de mélange plus proche de la réalité physique du mélange, que la matrice de mélange obtenue avec les méthodes comprenant une étape de blanchiment. Des aspects préférés mais non limitatifs du procédé selon l'invention sont les suivants : - les matrices de rotation de Givens comportent des éléments correspondant à des fonctions trigonométriques de 8 et les matrices de rotation hyperbolique comportent des éléments correspondant à des fonctions hyperboliques de cp ; - l'étape de diagonalisation comprend une étape de multiplication d'une famille de matrices estimées (' k) par le produit d'une matrice de rotation de Givens et d'une matrice de rotation hyperbolique à droite et la transposée de ce produit à gauche, et à calculer les valeurs de 8n et de cpn pour minimiser les termes non diagonaux des matrices résultantes d'ordre n de cette multiplication ; - l'étape de multiplication est réitérée sur les matrices résultantes d'ordre n ; - l'étape de multiplication est réitérée jusqu'à ce que l'angle de Givens 8 et � l'angle hyperbolique � (p soient tels que cose et coshcp tendent respectivement vers 1; - les matrices estimées Rk sont (à l'erreur d'estimation près) égales au produit de la matrice de mélange par une matrice diagonale par la transposée de la matrice de mélange (Rk = ADkAT ), ce sont par exemple des matrices de cumulants, ou des matrices de covariance dans différentes bandes de fréquence, ou des matrices d'intercorrélation avec différents décalages temporels, ou des matrices de covariance estimées sur différents intervalles de temps ; - les matrices de rotation de Givens comprennent des éléments en cose sur la diagonale, un élément en sine sous la diagonale et un élément en -sine sur la diagonale, et les matrices de � rotation hyperboliques � comprennent des éléments en coshcp sur la diagonale, un élément en sinh(p sous la diagonale, et un élément en sinh(p sur la diagonale ; 'cos° ùsine et les sine cose i "cosh(p sinh (p sinh (p cosh (p ~ L'invention concerne également un dispositif de séparation de signaux mélangés enregistrés par au moins deux capteurs en une pluralité de signaux composants, chaque signal mélangé correspondant à une combinaison linéaire instantanée des signaux composants issus de sources indépendantes, chaque signal mélangé étant apte à s'écrire comme le produit d'une matrice de mélange multipliée par un vecteur des signaux composants issus des sources, caractérisé en ce que le dispositif comprend en outre des moyens pour la diagonalisation conjointe dans lesquels la matrice de mélange est estimée à partir d'un produit de matrices de rotation correspondant à une ou des matrices de rotation de Givens et à une ou des matrices de rotation hyperbolique. L'invention concerne également un produit programme d'ordinateur comprenant des instructions de code programme enregistré sur un support utilisable dans un ordinateur, caractérisé en ce qu'il comprend des instructions pour la mise en oeuvre du procédé décrit ci-dessus. les matrices de rotation de Givens sont du type matrices de � rotation hyperbolique � sont du type PRESENTATION DES FIGURES D'autres caractéristiques, buts et avantages de la présente invention ressortiront encore de la description qui suit, laquelle est purement illustrative et non limitative et doit être lue en regard de la figure 1 annexée qui illustre un mode de réalisation du procédé selon l'invention. DESCRIPTION DE L'INVENTION On va maintenant décrire plus en détail un mode de réalisation du procédé selon l'invention en référence à la figure 1. Décomposer un mélange en composantes indépendantes est intéressant car il facilite l'interprétation, l'analyse ou le traitement de chacun des phénomènes qui se masquaient les uns les autres dans le mélange. Il peut s'agir par exemple : - des divers phénomènes électriques (cérébraux et musculaires) contribuant à un signal électro-encéphalographique (EEG), - des différents émetteurs radio dans une bande étroite de fréquence, - des émissions lumineuses de différents fluorochromes présents dans une solution, - des spectres de masses de différents mélanges des mêmes molécules, La séparation des différentes contributions permet alors : - d'étudier le signal EEG débarrassé des parasites musculaires (clignements, ECG, etc.), - de démoduler séparément le signal de chaque source radiofréquence, - d'estimer la concentration de chaque fluorochrome (cytométrie par exemple) ou 6 25 30 - d'identifier les différentes molécules présentes dans le milieu étudié par spectrométrie de masse. Le signal d'une source est fonction d'une grandeur variable qui peut être le temps (EEG, radiofréquence), la longueur d'onde (fluorescence), la masse (spectrométrie de masse), etc. Dans la suite, on présentera le procédé selon l'invention en considérant que les signaux des sources sont fonction du temps sans que cela restreigne la généralité du procédé. Le problème de séparation de sources se modélise mathématiquement de la façon suivante. Chacune des M sources émet un signal sm(t) indépendant des signaux émis par les autres sources. Le mélange (à un instant donné t) de ces M contributions est observé par un réseau de N capteurs (où N>M) qui fournissent N mélanges distincts. Les données observées à un instant t, notées x(t), sont donc multidimensionnelles. Le vecteur de ces N observations à un instant donné peut s'écrire comme le produit d'une matrice A (N lignes, M colonnes) dite � matrice de mélange � par le vecteur s(t) des signaux sources, soit : ^ a1M ( S1 (t) a2M s2(t) ~xN(t)/ ~aNl aN2 .. ^ aNM \sM (t) soit x(t) = As(t) xl (t) ail a12 x2 (t) a21 a22 Les contributions sm(t) des différentes sources aux mesures x(t) peuvent être mises en évidence en reformulant l'équation précédente : x1 (t) alm x2 (t) m=M a2m m=M Sm (t) = L amsm (t) m=1 \,xN(t)/ \aNm Le procédé selon l'invention permet d'estimer cette matrice de mélange A sans information a priori sur cette matrice (on la suppose seulement non singulière) puis à estimer les signaux sources. Les signaux sources sont supposés indépendants et on exploite certaines de leurs propriétés statistiques comme leur caractère nongaussien, leurs différences spectrales ou leur caractère non-stationnaire par exemple. Dans une première étape 100, on calcule une famille de K matrices estimées Rk (pour k compris entre 1 et K) à partir des données enregistrées par le réseau de capteur. Ces matrices estimées Rk sont égales (à l'erreur d'estimation près) à des matrices Rk qui ont la structure suivante : m=M Rk = ADk AT = L Dk (m, m) aman, m=1 où les Dk sont des matrices diagonales, et où AT est la transposée de la matrice de mélange A. Bien entendu l'homme du métier comprendra que les matrices Dk peuvent être diagonales par bloc. Le 3ième terme de l'équation n'est qu'une reformulation du 2nd qui fait apparaître les matrices Rk comme des sommes pondérées de M matrices de rang 1 (lesamam ). Le terme Dk(m,m)amam est la contribution de la source m à la matriceRk . Selon les contextes, ces matrices estimées Rk peuvent être des matrices de cumulants d'ordre 3 ou 4 (méthode JADE [1], [lb]), des matrices de covariance dans différentes bandes de fréquence [3], des matrices d'inter- =L m=1 corrélation avec différents décalages temporels (méthode SOBI [2]) ou des matrices de covariance estimées sur différents intervalles de temps [3]. Bien entendu, ces matrices estimées peuvent être d'autres types connus de l'homme du métier. Dans une deuxième étape 200, on effectue une diagonalisation conjointe des matrices estimées, comme décrit plus en détail dans la suite. Le procédé selon l'invention permet de calculer la matrice A à partir des matrices RI ,...,Rk,...,RK . On cherche une matrice B telles que toutes les matrices BTRkB soient � les plus diagonales � possible : BTRkB=Dk On construit B comme produit de matrices NxN élémentaires de la forme suivante : M(i, j,e,(p) _ a R II-i-1 j Y 8 où i 0 et ùx2+y'+z2 =1, l (i, i) ù 1 (j, j))/2 (i, j) (R~(i,i)ùRk(j,j))/2 (i, j) (RK(i,i)-RK(j,j))/2 RK(i,j), x 2 - calculer coshcp =l+ ~l+x2 et sinhcp = ~~1+.~l+x2 2cosh(p x - calculer cos0 = j z . j1+x2 2 1+ et sin 0 = ùy ( 2(1+x2) 1+ z .I1+x2 L'homme du métier appréciera que d'autres formules équivalentes existent puisqu'il ne s'agit que de résoudre les équations suivantes : 2 coshcp sinh cp = x ù 2 sin 0 cos0 (cosh2 p + sinh 2 (p) = y (2cos20 -1)(cosh2 p+sinh2 p) = z Les formules ci-dessus ont l'avantage d'avoir une bonne stabilité numérique quand x tend vers 0 en fin de convergence du procédé. D'autres modifications mineures peuvent être utiles pour s'adapter à la précision finie des calculs informatiques. Par exemple, la plus petite valeur propre positive théoriquement peut être, dans la pratique, très légèrement négative ; c'est néanmoins elle qu'il convient de retenir. Ces points sont évidents pour l'homme du métier. Le procédé comprend donc l'itération des produits : M(i,>,e,(p)T XRk xM(i, pour i La matrice de mélange A est alors égale à l'inverse de ML ((L)xGr (8L)x...xHo ((o)xGô (80) . Une fois la matrice de mélange A déterminée, il est possible de séparer les signaux sources. ùy 2 cos0 V1+ x2 Un exemple d'implantation en langage Matlab est donné ci-après : function [A,number_of_loops,matrices] = J_Di(matrices,threshold) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% % % J Di.m written the 19th of december 2007 by Antoine Souloumiac. % % This software computes the Joint Diagonalisation of a set of % matrices. % % input : % matrices: array of matrices to be jointly diagoonalized % threshold: threshold for stop criterium on maximal sines % % outputs: % A : estimated mixture matrix % number_of_loops: number of loops necessary for convergence % matrices: array of matrices after joint diagonalization % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% % initialization mat_size = size(matrices,l); nb mat = size(matrices,3); A = eye(mat_size); go_on = 1; number_of_loops = -1; % joint dyadization loops while go_on for n=2:mat_size for m=l:n-1 % initialization s_max = 0.0; sh max = 0.0; % Givens and hyperbolic angles calculation C = [ (squeeze(matrices(m,m,:)) + squeeze(matrices(n,n,:)))/2 (squeeze(matrices(m,m,:)) - squeeze(matrices(n,n,:)))/2 squeeze(matrices(m,n,:)) ]; [V,D] = eig(C'*C,diag([-1 1 1])); % normalisation (-1 1 1) of generalized eigenvectors V = V*diag(l./sgrt(abs(diag(V'*diag([-1 1 1])*V)))); % computation of the medium generalized eigenvalue [sortD,index] = sort(diag(D)); ind_min = index(2); v = V(:,ind min)*sign(V(3,ind min)); % computation of trigonometric and hyperbolic cosines and sines ch = sgrt( (l+sgrt(l+v(1)^2))/2 ); sh = v(l)/( 2*ch ) ; c = sgrt( ( 1 + v(3)/sgrt(l+v(1)^2) )/2 ); s = -v(2)/( 2*c*sgrt( 1+v(l)A2 ) ); % maximal sine and sinh computation for stop test s_max = max(s_max,abs(s)); sh max = max(sh max,abs(sh)); % product of Givens and hyperbolic rotations rmll = c*ch - s*sh; 35 40 45 50 55 60 rm12 = c*sh - s*ch; rm21 = c*sh + s*ch; rm22 = c*ch + s*sh; % matrices Givens update h_slicel = squeeze(matrices(m,:,:)); h_slice2 = squeeze(matrices(n,:,:)); matrices(m,:,:) = rmll*h_slicel + rm2l*h_slice2; matrices(n,:,:) = rm12*h_slicel + rm22*h_slice2; v_slicel = squeeze(matrices(:,m,:)); v_slice2 = squeeze(matrices(:,n,:)); matrices(:,m,:) = rmil*v slicel + rm2l*v_slice2; matrices(:,n,:) = rml2*v slicel + rm22*v slice2; % dyads Givens update coll = A(:,m); col2 = A(:,n); A(:,m) = rm22*colt -rm12*col2; A(:,n) = -rm2l*coll + rmll*col2; end end number of_loops=number of_loops+l; % renormalization of the A columns (equal norm) col_norm = sgrt(sum(A.^2,1)); d = (prod(col_norm).^(1/matsize))./col_norm; A = A.*repmat(d,mat_size,l)_ matrices = matrices.*repmat((1./d)'*(1./d),[1 1 nb mat]); % stop test if smax 100 go_on = 0; warning('maximum number of loops has been reached') end 40 end Le lecteur appréciera que de nombreuses modifications peuvent être apportées à l'invention telle que décrite ci-dessus sans sortir matériellement 45 des enseignements du présent document. Par exemple le produit des matrices de rotation de Givens (G(9)) et hyperbolique (H(q))) peut correspondre soit à la multiplication d'une matrice de rotation de Givens par une matrice de rotation hyperbolique, soit au produit d'une matrice de rotation hyperbolique par une matrice de rotation de 50 Givens. 15 10 15 20 25 30 Par conséquent, toutes les modifications de ce type sont destinées à être incorporées à l'intérieur de la portée des revendications jointes. REFERENCES: [1] Comparaison de méthodes de séparations de sources, Souloumiac A., Cardoso J.-F., Proceedings of GRETSI 1991, Juan les Pins, Septembre 1991. [1 b] Blind beamforming for non Gaussian signais, Jean-François Cardoso 10 and Antoine Souloumiac, ln IEE Proceedings-F, 140(6):362-370, December 1993. [2] A blind source separation technique using second-order statistics, Belouchrani A., AbedMeraim K., Cardoso J.-F., Moulines E., IEEE TRANSACTIONS ON SIGNAL PROCESSING 45 (2): 434-444, FEB 1997. 15 [3] Joint Approximate Diagonalization of Positive Definite Hermitian Matrices, Pham D. T., SIAM J. on Matrix Anal. and Appl., 22(4): 1136-1152, 2001. [4] A Fast Algorithm for Joint Diagonalization with Non-orthogonal Transformations and its Application to Blind Source Separation, Ziehe 20 A., Laskov P., Nolte G., Müller K.-R., Journal of Machine Learning Research 5 (2004) 777-800. [5] Quadratic Optimization for Simultaneous Matrix Diagonalization, IEEE Trans. on Signal Processing, vol. 54, no. 9, pp. 3270-3278, September 2006. 25 [6] Nonorthogonal Joint Diagonalization Algorithm Based on Trigonometric Parametrization, Wang F., Liu Z., Zhang J., IEEE Trans. on Signal Processing, vol. 55, no. 11, pp. 5299-5308, November 2007. REVENDICATIONS 1 Procédé de séparation de signaux mélangés (xi(t), ..., xN(t)) enregistrés par au moins deux capteurs en une pluralité de signaux composants (si(t), ..., sM(t)), chaque signal mélangé xn(t) correspondant à une combinaison linéaire des signaux composants (si(t), ..., sM(t)) issus de sources indépendantes, chaque signal mélangé étant apte à s'écrire comme le produit d'une matrice de mélange (A) multiplié par un vecteur des signaux composants issus des sources (si(t), ..., sM(t)), caractérisé en ce que le procédé comprend en outre une étape de diagonalisation conjointe (300) dans laquelle la matrice de mélange (A) est estimée à partir d'un produit de matrices de rotation correspondant à une ou des matrices de rotation de Givens (G(6)) et à une ou des matrices de rotation hyperbolique (H(cp)). 2. Procédé selon la revendication 1, caractérisé en ce que les matrices de rotation de Givens comportent des éléments correspondant à des fonctions trigonométriques de e et les matrices de rotation hyperbolique comportent des éléments correspondant à des fonctions hyperboliques de (p. 3. Procédé selon la revendication 2, caractérisé en ce que l'étape de diagonalisation comprend la multiplication d'une famille de matrices 25 estimées (Rk) par le produit d'une matrice de rotation de Givens et d'une matrice de rotation hyperbolique d'une part et la transposée de ce produit d'autre part, et l'estimation des valeurs de en et de cpn permettant de minimiser les termes non diagonaux des matrices résultantes d'ordre n de cette multiplication. 10 15 20 30 54. 5. 6. 15 7. 20 25 8.Procédé selon la revendication 3, caractérisé en ce que l'étape de multiplication est réitérée sur les matrices résultantes d'ordre n. Procédé selon la revendication 4, caractérisé en ce que l'étape de multiplication est réitérée jusqu'à ce que l'angle de Givens en et l'angle hyperbolique (Pn soient tels que cosen et coshcpn tendent respectivement vers un ou les sinus correspondants vers zéro. Procédé selon l'une des revendications précédentes, caractérisé en ce que les matrices estimées sont égales au produit de la matrice de mélange par une matrice diagonale quelconque par la transposée de la matrice de mélange, ce sont par exemple des matrices de cumulants, ou des matrices de covariance dans différentes bandes de fréquence, ou des matrices d'intercorrélation avec différents décalages temporels, ou des matrices de covariance estimées sur différents intervalles de temps. Procédé selon l'une des revendications précédentes, caractérisé en ce que les matrices de rotation de Givens comprennent des termes en cose sur la diagonale, un terme en sine sous la diagonale et un terme en -sine sur la diagonale, et les matrices de rotation hyperboliques comprennent des éléments en coshcp sur la diagonale, un élément en sinhcp sous la diagonale, et un élément en sinhcp sur la diagonale. Procédé selon l'une des revendications précédentes, caractérisé en ce que les matrices de rotation de Givens sont du type "cose -sine et les matrices de rotation hyperboliques sont du type sin 0 cos0 i "cosh(p sinh (p sinh (p cosh (p ~ 10 159. Dispositif de séparation de signaux mélangés (xi(t), ..., xN(t)) enregistrés par au moins deux capteurs en une pluralité de signaux composants (si(t), sM(t)), chaque signal mélangé xn(t) correspondant à une combinaison linéaire des signaux composants (si(t), ..., sM(t)) issus de sources indépendantes, chaque signal mélangé étant apte à s'écrire comme le produit d'une matrice de mélange (A) multiplié par un vecteur des signaux composant issus des sources, caractérisé en ce que le dispositif comprend en outre des moyens pour la diagonalisation conjointe dans lesquels la matrice de mélange A est estimée à partir d'un produit de matrices de rotation correspondant à une ou des matrices de rotation de Givens et à une ou des matrices de rotation hyperbolique. 10. Produit programme d'ordinateur comprenant des instructions de code programme enregistré sur un support utilisable dans un ordinateur, caractérisé en ce qu'il comprend des instructions pour la mise en oeuvre du procédé selon l'une des revendications 1 à 8.