FR2644915A1 - Procede et dispositif d'analyse spectrale en temps reel de signaux instationnaires complexes - Google Patents

Abstract

Le procédé et le dispositif permettent notamment l'analyse spectrale de signaux fournis par les vélocimètres Doppler médicaux. On modélise des tranches successives du signal échantillonné, par un processus autorégressif d'ordre élevé, proche ou égal du maximum autorisé par la durée des tranches et la fréquence d'échantillonnage, en faisant une estimation des paramètres caractérisant le modèle par filtrage de Kalman rapide utilisant un algorithme qui est simplifié dans toute la limite autorisée par la stationnarité locale du modèle dans la tranche et on détermine la densité spectrale de puissance g f à partir des paramètres estimés.

Description

Procédé et dispositif d'analyse spectrale en temps réel de signaux

instationnaires complexes

L'invention concerne l'analyse spectrale adapta-

tive en temps réel de signaux instationnaires à valeurs complexes représentant un phénomène physique. L'invention est utilisable chaque fois qu'il est

nécessaire d'analyser en temps réel de tels signaux re-

présentables par un modèle paramétrique auto-régressif, en faisant l'hypothèse que les coefficients du modèle restent stationnaires sur des intervalles de temps courts. Elle s'applique notamment à tous les domaines de l'imagerie médicale, par exemple par RMN, et plus généralement chaque fois qu'une analyse spectrale doit être effectuée sur des signaux organiques, par exemple en électrocardiographie et électroencéphalographie. Elle trouve toutefois une application particulièrement importante dans l'analyse spectrale des signaux fournis par les vélocimétres Doppler ultra-sonores servant à étudier l'écoulement sanguin dans les cavités cardiaques

ou les vaisseaux.

A l'heure actuelle, l'analyse spectrale paramé-

trique est relativement peu utilisée en vélocimétrie Doppler et en échographie ultra-sonore médicale. Elle reste au stade du laboratoire. Elle consiste à modéliser

les signaux à l'aide d'un modèle paramétrique, généra-

lement du type auto-régressif à moyenne ajustée, dont les paramètres sont ajustés par optimisation à l'aide d'algorithmes du type moindres carrés. Cette approche a l'avantage de permettre de travailler sur des tranches

courtes du signal mais elle impose d'utiliser des mo-

dèles d'ordre faible et a en conséquence une résolution en fréquence limitée, ce qui est particulièrement gênant

pour l'analyse spectrale de signaux fortement instation-

naires. Par ailleurs, les méthodes connues ne traitent pas les signaux démodulés à valeurs complexes comme tels, mais traitent séparément chacune des composantes,

respectivement en phase et en quadrature.

Les appareils du commerce utilisent des méthodes d'analyse spectrale nonparamétriques du type analyse de Fourier pour pouvoir en effectuer la mise en oeuvre à l'aide d'algorithmes de transformation de Fourier rapide. Mais ces méthodes nécessitent de travailler sur

des tranches longues du signal pour obtenir une estima-

tion convenable de la représentation spatiale recherchée, ce qui est incompatible avec la nature des

instationnarités rencontrées dans les domaines précé-

demment cités.

L'invention vise notamment & fournir un procédé et un dispositif d'analyse spectrale répondant mieux que ceux antérieurement connus aux exigences de la pratique, notamment en ce qu'ils améliorent la résolution en fréquence, permettent une analyse plus précise des signaux et restreignent la variabilité des spectres obtenus pour représenter des écoulements dans des

conditions voisines.

Dans ce but, l'invention propose notamment un procédé d'analyse spectrale adaptative par modélisation du phénomène représenté par le signal, caractérisé en ce qu'on modélise des tranches successives du signal échantillonné, par un processus auto-régressif d'ordre élevé, proche du ou égal au maximum autorisé par la durée des tranches et la fréquence d'échantillonnage, en faisant une estimation des paramètres caractérisant le modèle par filtrage de Kalman rapide utilisant un algorithme qui est simplifié dans toute la limite autorisée par la stationnarité locale du modèle dans la tranche et en ce qu'on détermine la densité spectrale de

puissance à partir des paramètres estimés.

Pour éviter les risques d'instabilité liés au choix d'un ordre élevé, l'estimation des paramètres du modèle auto-régressif est régularisée par mise en oeuvre d'une contrainte consistant à minimiser un critère de la forme: J2(a) - Jl(a) + 0(a) o Jl(a) est la somme des moindres carrés [y(n)-z(n)]2, z désignant la fonction échantillonnée modélisant le signal par l'intermédiaire des paramètres a à déterminer et y la fonction représentée par les échantillons du signal et o 0(a) est une fonctionnelle régularisante tenant compte de la connaissance a priori de la nature

de la fonction ou obtenue par optimisation.

Grâce à l'utilisation d'un modèle d'ordre élevé, on obtient une très bonne définition en fréquence. Le procédé permet d'obtenir directement des paramètres de calcul de la densité spectrale de puissance du signal

dans la fenêtre d'analyse.

L'invention propose également un dispositif d'analyse spectrale directement adaptable sur un vélocimètre Doppler classique, sans modification de la sonde existante et de son électronique associée, ni des

moyens de visualisation.

L'invention sera mieux comprise à la lecture de

la description qui suit du procédé et du dispositif

selon l'invention et du cheminement qui a permis d'y aboutir. Elle se réfère aux dessins qui l'accompagnent, dans lesquels: - la' Figure 1 est un schéma de principe montrant les paramètres qui interviennent en vélocimétrie Doppler du flux sanguin dans le cas particulier d'une technique de Doppler continu; - la Figure 2 est un schéma montrant l'enveloppe d'un spectre Doppler représentatif; - la Figure 3 est un schéma de principe de la démodulation du signal Doppler; - les Figures 4 et 5 sont des synoptiques montrant une implémentation possible du procédé suivant l'invention. Modélisation des signaux non-stationnaires

Avant d'aborder la description de l'invention,

il convient de rappeler des indications sur la modélisa-

tion des signaux représentant des phénomènes non sta-

tionnaires et sur la régularisation du modèle.

Dans le cas o la non-stationnarité du signal est suffisamment lente devant la vitesse de variation de l'amplitude du signal, celui-ci peut être considéré

comme "localement stationnaire", c'est-à-dire station-

naire sur une tranche d'analyse de N échantillons, N étant typiquement de l'ordre de quelques dizaines. Les échantillons observés y(n), n t 1, 2,... , N, sont décrits par un modèle auto-régressif (AR) exprimé par: y(n) = z(n)+b(n) (1) p z(n) = Z aiy(n-i) = atyp(n) (2) i=1 avec at = [apapl,...., al] (3) yp(n)t = [y(n-p),..., y(n-1)] o at est le vecteur des coefficients du modèle dont l'ordre est p. et o b(n) représente le processus

générateur du signal, de variante a2 inconnue.

L'analyse spectrale s'effectue alors suivant un critère qui est souvent celui des moindres carrés, consistant à minimiser Jl(a): N Jl(a) = Z [y(n) - z(n)]2 (4) n-1 Le modèle étant linéaire, les algorithmes de calcul sont rapides. Mais l'obtention d'une bonne approximation est subordonnée au choix d'un ordre p aussi élevé que possible; or, si on prend p proche de N ou égal à lui, la minimisation de J1 n'est pas un critère satisfaisant par suite de la variance excessive

de la solution.

Pour écarter cette difficulté et stabiliser la solution, on a proposé de régulariser le modèle en ajoutant à J1 une fonctionnelle stabilisatrice, ce qui permet de conserver une valeur élevée de p (proche de ou même égale à N). Le critère consiste alors à minimiser J2(a): J2(a) J1(a) + 0(a) (5) o O(a) est une fonctionnelle régularisante, dont il

existe de nombreux exemples.

Le choix de la fonctionnelle régularisante n(a) traduit l'information a priori de la nature des propriétés spectrales locales du signal analysé. Cette information est que la densité spectrale de puissance

inconnue est une fonction "douce", c'est-à-dire pré-

sentant une certaine régularité, comparée à un spectre inorganisé constitué d'une multitude de raies sans relation entre elles. Pour des raisons de commodité de

mise en oeuvre, on choisit une fonctionnelle quadra-

tique, ce qui présente l'avantage de conduire à une solution explicite et linéaire par rapport aux valeurs

observées y(n).

Une mesure de "douceur spectrale" est constituée par 1/2

Dk = I 8akH(f) 12 df,k 5 0,1,2,...

-1/2 afk o H(f) est la fonction de transfert en fréquence du filtre blanchisseur du signal analysé, qui est aussi l'inverse de son filtre générateur qui est décrit par les équations (1) et (2). Une valeur élevée de Dk signifie une densité spectrale de puissance non douce, et ce d'autant plus que l'ordre k de la dérivée est élevé. On utilise une douceur spectrale d'ordre zéro

J1/2 P

Do I H(f) 12 dfa 1 + Z a2 1/2 m=l qui est, & une constante près, la norme euclidienne du vecteur des paramètres a, ce qui conduit à adopter la fonctionnelle régularisante P 2 0(a) = Z a.. + (6) m ma=l o le coefficient p est un nombre réel positif appelé "coefficient de régularisation" et dont il sera vu qu'il Joue un très grand rôle dans. le fonctionnement de l'invention.

Modélisation & régularisation suivant l'invention.

L'invention utilise une démarche qui peut étre considérée comme en deux étapes, utilisant notamment les constatations que l'on a très souvent affaire, dans la pratique, & des signaux pour lesquels on peut effectuer une étude locale du phénomène par division de l'horizon temporel en blocs de longueur N, et qu'on peut faire l'hypothèse d'une stationnarité locale dans chaque bloc,

de manière à permettre d'exploiter des propriétés d'in-

variance locale du modèle, en effectuant donc une analyse spectrale adaptative reflétant une variation lente des paramètres a vis-à-vis des fluctuations d'amplitude du signal, à savoir: 1. On recherche un modèle localement stationnaire, c'est-A-dire un modèle invariant équivalent défini par

l'optimum du critère (5) calculé sur une fenêtre tem-

porelle réduite, en effectuant un compromis entre une trop grande longueur de fenêtre (qui ne permet pas de suivre fidèlement les non- stationnarités du signal) et une trop courte longueur (qui ne permet pas d'estimer le modèle dans de bonnes conditions), ce qui conduit à

choisir un modèle AR de l'ordre p le plus élevé pos-

sible, proche de N ou égal à N. imposant de stabiliser

le modèle.

Ce modèle est recherché pour la valeur du coefficient de régularisation m et celle de la variance c2 du processus générateur qui rendent maximale la

vraisemblance du modèle.

2. On simplifie les calculs de maximisation de vraisemblance du modèle en adoptant une solution sous-optimale: - on réduit le domaine de variation de w de l'équation (6) du critère de minimisation à un ensemble fini et présélectionné de valeurs discrètes, et - on calcule, par un algorithme qui peut être rapide,

la solution optimale à et éventuellement la vrai-

semblance correspondante V(cr2, s/y) à partir des échantillons du signal observé, pour retenir la plus élevée. s La première étape conduit à organiser les échantillons du signal en blocs consécutifs de longueur N et & effectuer sur chaque bloc une régularisation. On impose une contrainte liant localement l'ensemble des paramètres et traduisant une information de douceur 1o spectrale: cela conduit à minimiser, sur chacun des blocs, un critère de la forme: J2(a)= Jl1(a) +. la-a0 12 (7) Lorsque le processus générateur du signal analysé est supposé gaussien, cela revient à admettre que les paramètres a sont caractérisés, dans une interprétation bayésienne, par la loi de probabilité a priori: f (a / ao, l,u2).Cexp là la-a0121 L'initialisation du vecteur des paramètres AR se

fait en choisissant ao = 0 dans le cas d'un signal sta-

tionnaire. On verra plus loin que cette façon de procé-

der peut être conservée pour certains signaux non sta-

tionnaires. Dans les autres cas on choisira, pour chaque

bloc courant, le résultat du traitement du bloc précé-

dent lorsque les blocs sont adjacents.

Il reste encore à fixer la valeur de p et de "2

pour pouvoir effectuer la minimisation du critère (7).

Ce choix est fait, comme il a été dit plus haut, en maximisant la vraisemblance du modèle. On verra que p et 42, appelés hyperparamétres, peuvent étre découplés dans

cette maximisation.

Pour exposer la première étape, il convient de donner les équations d'état associées au problème et de montrer comment on peut substituer au filtre de Kalman standard (qui fait appel & une équation de Riccati n'exploitant pas la propriété de stationnarité locale du modèle qui est à la base des méthodes adaptatives, ni la propriété de décalage des vecteurs yp(n) successifs qui est la conséquence du choix d'un modèle AR) un filtre plus simple, mettant en oeuvre les équations de Chandrasekhar, pour le calcul des paramètres à qui minimisent, p et o2 étant fixés, le critère (7), ainsi

que pour le calcul de la vraisemblance correspondante.

Pour cela, on fait apparaître explicitement ces propriétés en définissant un "vecteur-paramètre étendu" ou vecteur d'état à m composantes: am(i) l [0ti- ap _l,.... a1,,O,,..., O]t, avec m>p et on considère le vecteur global des données: Ym = [y(-p+l),y(-p+2),..., y(O), y(1),..., y(m-p)lt Le modèle linéaire défini par les équations (1) et (2) peut alors se mettre sous la forme du modèle d'état équivalent suivant:

aó(i+l) - Dam(i) i. 1,2,...

y(i) - ytm am(i)+b(i) o D est l'opérateur de décalage:

O.......

1 0....

D 0 i 0....

0 0 1 0

Le calcul récurrent de la solution peut alors être effectué par le filtre de Kalman suivant: &m( i+l/i)=D&m( i/i-1)+( i)r i)-1[y( i)-ytmam( i/i-)] ( 8 r(i) - ytmP( i/i1-)ym + a2 km(i) = D Pm(i/i-I)ymi Pm(i+1/i) = D Pm(i/i-l) Dt - km(i)r( i)-lk*m(i) avec les conditions initiales: &m(1/O) = [a(1/0), o..., 0]t r(l) = ytm(1/O)ym + 2 (10) KM(l) = D Pm(l/O)Ym P(1/o) 2 Pm(1/0) =......._ p(1/O)=,2 r= " 0!0 ou Ym désigne le vecteur simplement conjugué du vecteur

Ym et k* désigne la matrice transposée conjuguée de k.

En apparence, l'introduction de ces variables supplémentaires est défavorable, car elle augmente les dimensions du vecteur d'état. Mais cette impression est

trompeuse, car cette introduction fait apparaître expli-

citement les propriétés d'invariance que l'on exploite.

On peut d'ailleurs vérifier qu'à chaque récur-

rence, le nombre de coordonnées non-nulles du vecteur gain k(i) reste égal à p. Ce modèle d'état étant invariant (D et Ym étant constant) et à variance stationnaire (a2 étant constant sur un bloc), on obtient directement un algorithme rapide en remplaçant les équations standard par des équations de Chandrasekhar. Non seulement ces équations permettent une résolution rapide, mais encore elles

permettent de passer, de façon naturelle, à des algo-

rithmes en racine carrée offrant une bonne stabilité numérique. On désignera par la suite par "Algorithme B-CAR"

les règles de régularisation mettant en oeuvre l'algo-

rithme,appliquées à des blocs de taille appropriée à la nature du signal. Ces équations ne propagent que les incréments des quantités nominales de l'équation (9), qui sont eux-mêmes factorisés. L'application directe de ces techniques de factorisation au filtre (8) (9) conduit aux équations suivantes: vm(i)t YmtB(i) wm(i)t - vm(i)t M(i) M(i+l) - M(i) + wm(i)r(i)lwm(i)t (11) km(i+l) - km(i) + D B(i)gm(i) r(i+l) - r(i) + vm(i)tm(i) B(i+ l) - D B(i)-km(i+l)r(i+l)-lvm(i)t On utilise dans ces équations le fait que M et Pm sont des matrices à symétrie hermitienne: M*=M et P*m=Pm. Une écriture condensée de l'algorithme est obtenue par projection, en ne faisant apparaître que les composantes intervenant effectivement à chaque récurrence dans le calcul. On définit pour cela un vecteur-gain effectif k(i), de dimension p. contenant les coordonnées non-nulles de km(i): km(i) - JOti+l, k(i)t, Otm-p-i-lIt De la méme façon, si a est la dimension de M(i), on peut aussi écrire: B(i) = 10a.,k, S(i)t, Oa,m-p-i-_lt

ou S(i) est une matrice de dimensions [(p+l),a. L'algo-

rithme peut alors se mettre sous la forme v(i)t = [y(i-p), yp(i+l)t]S(i) w(i)t = v(i)t M(i) M(i+l) = M(i) + w(i)r(l)'lw(i)t r(i+l) = r(i)+v(i)t w(i).(12) O k<i) à..=L.... + s(i)(i) k(i+1 S(i+l) S(i) - |-..---- ri+l)lv(i)t k( i+)p la remise à jour des paramètres se faisant selon: [t.. 'o, ',s....Z-dA,'dO]," p

'O'O'T'O ''],.A 0C

[.... 'o'o 'IT-dA. ZA '0]." 6 o'o'r] = m-s 90 (UA twp + p + UA)OA + ( 's + U Ms + ms 'Ws)OA-. (O/T);d-a (O/1)mdaSZ Jssoduoonp ind uo sOlQe & -(O/T)d

A I A

*-- ZA 1A OZ

asod uo, Ts 'Ieigug6 seo aT Suva 0! 0].(o)m 0O (0/1)d (O/1)md -(T)ux.T) T) r(T)u!-{a (O/T)mda - (Z)md9 : enpue, aw.xo. snos suoTpenb9 xne ueueaasz ua 'lzTzOg.s zuewalTe6 D.nad aoueTzsaoo ap eT:TUT,usmaXouT" (I)g(t)H()(O/l)"ld-(I:Z)Nd- (Z)Md : TeT:.TUTuemouT,T 0 ap UoTs.esTaozoe; el ap s0ol (I)w 4.S (I)B 4uUsTnpgp es wuop (O/1) d TboTxd e eoueTaeAoo sp SoTz;nm el zed sau

-UOTZTPUOD zsa am4qnT206teT Sp UOTFleSTfITZTUT.

(eXTeUUOTn.e.s-uou seo el suep susosepe sooIq Sp.uem -azTe'a) aoueTZeAoo,; <dAz np seuzgTqoad sep,nb (<STOU ç -uoTzel.s-uou seo Sl suep 9TosT ooTq un,p no SaTmUuOTZess seo np z.ueamreaza) sazuse;a.xd. sswgTqoa.d ssp usTq Tssne zanbTTdde.s sp Se uzuAU.T e SwuqTaOBTQs ze (ET) [(T-I).(T)d- () ,]_(A)A(T)+()T-z/)-(/ T+T)> St La matrice P(1/O) est choisie diagonale, et le rang de amP(2) est en général "=3. La factorisation initiale n'étant pas unique, on choisit avantageusement pour M(1) la matrice de signature de &Pm(2) qui fournit l'amorce d'un algorithme en racine carrée dont les

qualités numériques sont meilleures.

Si la matrice de covariance a priori P(1/O) est diagonale, les vecteurs gm et dm sont identiquement nuls et l'incrément initial &Pm(2) se réduit à: 6Pm(2) = -vOSm S*m + vo Vm v m - km(1)r(1)-lkm(1)*, de rang au plus égal à 3 puisqu'il se décompose en une somme de trois matrices dyadiques ou antiscalaires, chacune de rang égal à 1. On peut écrire &Pm(2) sous la forme: 6Pm(2) = B(1)M(1)B(1)* avec comme facteur central:

-1 0 O]

M(1) = -1

qui conduit, une fois la projection effectuée pour obte-

nir l'algorithme en dimension réduite, à: S(1) = [k(1),b,c] k(l) = [O,k(l) tr(1)-'/2']t k(l) = yp(1) b = c = [0,0,0,...,,13t r(1> ytmPm(l/O)Ym + o- = yptl)typ(1)+G' sous les conditions, généralement remplies, que: - les coefficients ai du modèle, qui sont complexes, aient des parties réelles et imaginaires dont les lois de probabilité a priori sont normales, centrées, indépendantes, et de variance (1/2)T2 (on retrouve ainsi la matrice de covariance a priori P(1/o)ur21), - le bruit générateur du signal observé b(n), qui est aussi complexe, ait des parties réelle et imaginaire dont les lois de probabilité a priori sont aussi normales, centrées, indépendantes, et de variance

(1/2)"2.

Cette solution rigoureuse utilise les valeurs effectivement observées antérieures à l'instant initial de la fenêtre dans-le vecteur général Ym. Cette méthode, Io que l'on peut qualifier de "covariance" suppose que l'on travaille: - soit sur un bloc de données isolé de longueur double de l'ordrep choisi, - soit sur un bloc de données adjacent à un bloc

antérieur dont le contenu a été conservé.

Le décompte des opérations arithmétiques élémen-

taires, limité aux seules multiplications, montre que la

complexité de cet algorithme est e(11p) par récurrence.

Si on choisit p-N, on voit que la complexité totale pour traiter N données est 9(N2), ce qui est normal pour un

algorithme de ce type.

Une solution permettant de simplifier l'algo-

rithme peut être qualifiée de "préfenétrage". Elle fait l'hypothèse que les p valeurs antérieures à l'instant initial dans la fenêtre de mesure sont nulles: avec pour conséquences: - <(l). D Pm(l/O)Ym ' O r(1) yp(1)t ip(1)+,' 6Pr(2) - yOSme m+VOVmv m Le rang de l'incrément initial est a=2. On peut choisir comme facteur initial: M(l) - [ 0 l 1 J [-O l J S(1) - t b,c]

On réduit ainsi sensiblement la complexité du calcul.

Lorsque la nature du phénomène représenté par le signal est bien connue, les hyperparamêtres peuvent étre choisis a priori. Dans le cas contraire, un calcul de la

vraisemblance des hyperparamntres est utile.

L'utilisation d'un filtre de Kalman permet de

calculer par récurrence la vraisemblance des hyperpara-

mètres (ou son logarithme) à l'aide des quantités

intervenant dans l'algorithme.

La recherche du maximum de vraisemblance revient en fait à un problème de moindres-carris régularisés déjà représentés par l'équation (7). Ce problème a une interprétation bayésienne simple puisque la minimisation de (7) est équivalente à la maximisation du critère de vraisemblance V(a):

{ 1 N

V(a) a exp t- (Z [y<n)-z(n)]2+ Z a2m (153 2e n=l m=l Lorsqu'on fait intervenir les hyperparamètres p. le calcul montre que la solution à est en fait la

moyenne de la loi a posteriori définie par la distribu-

tion conditionnelle des données:

2 { 1 N23

f(y/a, ')o.exp - (n [y(n)-z(n)]2) n=! et la loi a priori: f(a/,,a) eexp {_ ata} 2a On peut donc considérer que cette loi a priori spécifie une classe d'estimateurs par l'intermédiaire des paramètres p et a'. Ce sont les hyperparamètres du problème. Comme les lois sont normales, la vraisemblance des hyperparamétres est représentée par une intégrale: V(o, ',p,y) - I f(y/a,u',)f(a/.pr')da (16) Le modèle à retenir est celui qui maximise cette vraisemblance par rapport aux hyperparamètres. En ce sens, la méthode est adaptative puisque le choix de valeurs a priori peut dépendre des données. Enfin, par application séquentielle de la règle de Bayes à la relation, on obtient: N V(", ""/y) - f[y(l)] n f[y(n)/y(l),..., y(n-l)] n=2 On utilise la densité marginale conditionnelie

de y(n) que permettent d'obtenir les échantillons obser-

vés y(l), y(2),...y(n-1): fty(n)/y(l),--y(n-l) -ifty(n)/a] fta/y(l),..., y(n-l)]da 1 1 p. : 2._ exp {[y(n) - i a(i)y(n-i) 2e Le calcul de cette densité conditionnelle exige S de connaître u et a. Mais y(n) - yp(n) t&(n/n-1)+e(n)

o &(n/n-l) est l'estimée optimale (c'est-a-dire condi-

tionnelle au passé (y(l), y(2),..., y(n-l)), et o e(n) est l'innovation du processus observé. On a donc aussi: e(n)2 f[y(n) / y(l),.., y(n-l)]ex.2.() Reste à évaluer la variance de l'innovation es(n)2. Mais un filtre de Kalman n'est sensible qu'à la

seule variance relative du bruit d'observation (c'est-

à-dire ici du processus générateur) et du bruit d'état

(c'est-à-dire ici de la variance a priori des para-

mètres). Le calcul récurrent de à(n/n-1) ne dépendant donc que de à et pas de la valeur de, on peut prendre I arbitrairement oe = - i dans l'algorithme; on a alors: re(n)2 = o2r(n) - aos[yt Pm(n/n-l)ym+l] d'o l'on tire: fy/y1.... f e(n)2}

f[y<n) / y<l),..., y(n-1)] = exp..

2 ( n3 t 2,2r(n) 0 Cette formule montre que les deux hyperparamètres sont A n

découplés. L'estimée: de 2 selon le maximum de vrai-

semblance (marginale) est donc: 1 N e(i)2

C= -z -

N i=l r(i) et la vraisemblance (marginale) de " est alors: 1 N e(i)2 L(p/y) = - f - + înlr(i)>] +cste 2 i=l r(i) En résumé le procédé complet comprend, lorsqu'on doit déterminer les hyperparamêtres optimaux: - pour chacune des valeurs discrètes de p choisies pour le calcul de J2, l'estimation d'un vecteur 'â(N/N) avec l'algorithme (10) (11), avec "* = 1; - le choix des paramètres & et des hyperparamètres "2 et y correspondant à la vraisemblance la plus forte; et - le calcul de la densité spectrale de puissance g(f) correspondante & l'aide de la relation classique p g(f) - es l-Z ak exp(-2ifk)12] (17) k. I pour - 1/2 À f C + 1/2 Application de la modélisation suivant l'invention aux signaux de vélocimètrie sanguine Doppler Le choix de la durée des fenêtres ou blocs de modèle auto-régressif, du nombre d'échantillons à prélever par bloc, de l'ordre du modèle et du critère de régularisation implique de connaître les caractères principaux des signaux Doppler. Ces signaux sont obtenus par un montage du genre schématisé en Figure 1. Le champ de vitesse étudié est éclairé par une onde acoustique monochromatique ou pulsée provenant de l'émetteur 10 et on mesure l'onde réfléchie par les cibles matérielles (hématies) se déplaçant sous l'action de ce champ de vitesse, recueillie par un récepteur 12 comportant un démodulateur. Si le champ de vitesse était stationnaire et uniforme, le signal reçu r(t) se déduirait du signal émis a(t) par une simple translation en fréquence. Ce n'est pas le cas, du fait de la contribution d'un grand nombre de cibles élémentaires présentant une dispersion par rapport & une cible idéale moyenne, le signal reçu a un spectre continu, de largeur de bande étroite par rapport à la fréquence de la porteuse (quelques kHz par rapport à quelques MHz). Dans le cas fréquent d'un signal Doppler à émission pulsée, la forme de l'onde émise entre également en Jeu. On obtient dans les deux cas un spectre & bande étroite, du genre montré en figure 2, se limitant à deux bandes symétriques de largeur 2F centrées en -fo et fo respectivement, F étant

très inférieur & la fréquence fo de la porteuse.

Pour simplifier le traitement des signaux, les appareils existants utilisent une démodulation, généralement synchrone, qui conduit à décomposer le signal reçu r(t) en deux composantes, respectivement en phase p(t) et en quadrature q(t). Ces signaux à basse fréquence (quelques kHz) n'exigent qu'une cadence d'échantillonnage réduite car leur spectre utile est réduit par un facteur 1000 environ par rapport à la porteuse. Le montage peut être celui schématisé en

figure 3.

En représentation complexe, on peut écrire: e(t) = Vr2 R. {elt)exp(iu(t) o wO0 est la pulsation de la porteuse et o é(t) désigne l'enveloppe complexe de modulation du signal émis, et r(t) i -2 p(t)cos(wot)+ r2q(t) sin(wot) o p(t) et q(t) sont des signaux & basse fréquence générés par le démodulateur 14 de la figure 3, constituant respectivement la partie réelle et la partie

imaginaire du "signal Doppler" à valeurs complexes y(t).

Le calcul montre qu'on peut obtenir le spectre du signal reçu r(t) par transformation de Fourier de la covariance du signal Doppler y(t). Par ailleurs p(t) et q(t) sont en général corrêlées, ce qui interdit de faire l'hypothèse - pourtant admise dans le traitement habituel du signal Doppler - que la covariance mutuelle Rp0(4g) de p et q ou leur densité interspectrale de puissance GpQ(f) est nulle. Or cette hypothèse est

nécessaire si on traite p(t) et q(t) séparément.

L'invention traite donc simultanément p(t) et q(t) ce qui est possible du fait que la méthode d'analyse spectrale décrite plus haut a été Justement conçue pour s'appliquer à des signaux à valeurs complexes. Choix des paramètres de réglage dans l'application à la vélocimétrie Doppler

Trois paramètres de réglage doivent être choi-

sis, soit a priori (ce qui sera le cas général étant donné que l'on connaît la nature du signal Doppler),

soit à l'issue d'une procédure d'optimisation.

Il s'agit de: - N, durée d'une fenêtre d'analyse ou longueur d'un bloc de données; - p. ordre du modèle AR adopté pour décrire le signal Doppler, a, coefficient de régularisation ou hyperparamétre.

Le mode d'initialisation doit aussi être choisi.

Il est défini par: - le choix des paramètres initiaux a0 au début de chaque bloc: a0 peut être ou bien nul, ou bien le résultat du traitement du bloc précédent; - le choix entre la version dite "covariance" et celle

dite upréfenétrée".

Le choix des trois paramètres N, p et A résulte d'un compromis. Augmenter N et p améliore la qualité du modèle, mais on est limité dans cette voie par la non-

stationnarité du signal qui pousse à réduire N et par la charge de calcul qui pousse à réduire p. On obtient de bons résultats avec une fréquence d'échantillonnage de l'ordre de quelques dizaines de kHz

et avec des fenêtres d'observation d'une longueur maxi-

male de N - 64 points, c'est-à-dire d'une durée maximale de 1 &' 2 ms (au lieu de 5 à 10 ms dans les appareils commercialisés actuellement). Comme on n'a pas besoin

d'une analyse spectrale du signal Doppler à des inter-

valles de temps aussi brefs, on peut organiser les données échantillonnées en blocs de N = 64 points au maximum et ne travailler que sur un bloc sur dix. Ceci o0 permet d'obtenir une densité spectrale de puissance toutes les 10 à 20 ms au plus, c'est-à-dire à raison

d'une cinquantaine par période cardiaque.

L'initialisation peut être faite par préfené-

trage ou par covariance, comme indiqué plus haut. Le préfenêtrage, actuellement utilisé, diminue le nombre de

multiplications à effectuer, mais nuit à la vraisem-

blance du modèle, surtout si p est proche de N. La méthode de covariance réduit l'ordre maximal du modèle à p = N/2 si on commence les récurrences à l'instant p+l (en cas de fenêtres isolées) et elle est un peu plus complexe. Le choix de l'hyperparamètre p est essentiel car il règle le degré de douceur de la solution et doit maximiser la vraisemblance du modèle du signal par rapport aux données: les valeurs des paramètres AR calculés ne dépendent en effet pas des variances a2 et T2, mais uniquement de leur rapport l. L'expérience montre que la courbe de variation de la vraisemblance en fonction de p varie très peu au voisinage du maximum et que les résultats de l'analyse sont peu affectés par des

variations d'un ordre de grandeur de ce paramètre.

L'expérience montre aussi que, pour un signal Doppler donné, les variations de la valeur optimale de l d'un bloc de données à un autre sont aussi d'un ordre de grandeur. La complexité de la méthode peut donc être réduite en substituant, à l'étape de maximisation de la vraisemblance du modèle qui doit être effectuée sur

chaque bloc de données, une étape préalable de classi-

fication des signaux en fonction des valeurs constantes & attribuer à M. Pour que ce soit possible, il est indispensable de normaliser M par rapport à la puissance du signal Doppler. En effet, le coefficient intervient dans le critère minimisé: J2(a) - Jl(a) + pla-a012 (7) o Jl(a) a la dimension de l'énergie du signal et les paramètres AR sont sans dimension. Le coefficient M a donc une dimension et varie comme le carré de l'échelle

dans laquelle est mesuré le signal Doppler.

Pour rendre M indépendant de la puissance du signal, il faut donc choisir: M - M.o=(N/p) o mo0 est le nouveau paramètre de réglage, cette fois sans dimension, et us est la variance du signal qui peut être fournie par le module d'acquisition du dispositif utilisé (24). Cette normalisation permet de travailler avec O0 constant tout au long d'un cycle cardiaque, indépendamment des fluctuations de puissance du signal Doppler. L'invention est susceptible de nombreuses

variantes de réalisation. En particulier, on peut rédui-

re la sensibilité du procédé à la précision limitée des calculs en propageant les racines carrées des matrices de covariance plutôt que les matrices elles-mêmes: une telle modification s'effectue de façon simple, comme

indiqué plus haut.

Constitution du dispositif Le dispositif utilisé pour la mise en oeuvre du procédé peut avoir la constitution générale montrée en figure 4, qui recherche un compromis entre la qualité des résultats et le volume des calculs, est adaptable sur des appareils existants et utilise les possibilités intrinsèques de processeurs de traitement numérique du

signal actuellement disponibles.

Le dispositif comporte deux processeurs couples

pour effectuer les calculs du spectre Doppler. Le pre-

mier, 20, calcule les paramètres d'un modèle du signal Doppler reçu, par l'algorithme B-CAR. Le second, 22, fournit, à partir de ce modèle, le spectre. Ce second processeur 22 commande une carte d'acquisition 24 qui acquiert les deux voies Doppler simultanément, à une fréquence d'échantillonnage qui est par exemple de 40 kHz pour l'étude de la circulation sanguine, choisie pour optimiser le compromis entre la durée des tranches analysées et la précision fréquentielle de l'analyse spectrale. Pour éviter une éventuelle dissymétrie dans l'acquisition des signaux Doppler, un seul convertisseur est utilisé. La carte peut être prévue pour détecter une

éventuelle saturation des signaux Doppler.

Le processeur 20 est programmé pour mettre à l'échelle l'hyperparamétre p de la méthode à partir de l'estimation de la puissance instantanée du signal Doppler fournie par le processeur 22, et calculer les paramètres AR du signal Doppler sur chaque bloc de

données transmis par le même module.

Le processeur 22 calcule le spectre du signal sur un bloc, & partir des paramètres AR fournis par le

processeur 20 sous le contrôle du module de visualisa-

tioh 26 qui fixe le nombre de points sur l'axe des fréquences et le nombre de niveaux de gris; il peut être prévu pour calculer une éventuelle composante continue dans la chalne d'acquisition. Il doit de plus évaluer la puissance instantanée du signal, utilisée par

le processeur 20 pour optimiser les résultats fournis.

L'architecture de l'appareil peut être celle

montrée en Figure 5. Le processeur 22 exécute l'algo-

rithme de calcul de la densité spectrale de puissance et gère les interfaces 28-30-32 d'entrée et de sortie de données. Le processeur 20 exécute l'algorithme B-CAR. Un mini-ordinateur hôte constitue une machine maître qui charge le programme dans les processeurs et

qui commande l'exécution des programmes. Une implanta-

tion en ROM du programme permettrait d'éviter la 1o nécessité d'un miniordinateur hôte. Les interfaces d'acquisition 30 et de visualisation 32 transfèrent les données, d'une part entre le dispositif d'acquisition et le processeur 22, d'autre part entre le dispositif de visualisation 26 et ce même processeur 22. Chaque

processeur est muni d'une mémoire propre de travail 36.

Il est nécessaire d'effectuer les calculs en format 32 bits avec virgule flottante. Des circuits disponibles sur le marché permettent de répondre à de tels besoins. Il s'agit d'une part des processeurs de traitement du signal composés d'une unité de calcul 32 bits en virgule flottante, d'une unité de commande, de mémoires et de dispositifs d'interface intégrés sur un seul circuit (par exemple DSP32 de ATT et,PD77230 de NEC), d'autre part, des unités arithmétiques 32 bits. On

a notamment implémenté le procédé à l'aide de proces-

seurs DSP32 dont les bus internes sont accessibles de l'extérieur du boîtier (ce qui permet d'augmenter la taille mémoire de 4 Koctets à 56 Koctets), et qui ont une interface parallèle permettant d'effectuer des accès-mémoire en DMA quelle que soit l'activité du processeur et une interface série & haut débit (8 Mbits/sec).

Claims (7)

REVENDICATIONS

1.Procédé d'analyse spectrale en temps réel de

signaux instationnaires à valeurs complexes, représen-

tant un phénomène physique, suivant lequel on modélise le phénomène représenté par le signal, par un processus auto-régressif, caractérisé en ce qu'on modélise des tranches successives du signal échantillonné, par un processus auto-régressif d'ordre élevé, proche ou égal du maximum autorisé par la durée des tranches et la o10 fréquence d'échantillonnage, en faisant une estimation des paramètres caractérisant le modèle par filtrage de Kalman rapide utilisant un algorithme qui est simplifié dans toute la limite autorisée par la stationnarité

locale du modèle dans la tranche et en ce qu'on détermi-

ne la densité spectrale de puissance g (f) à partir des

paramètres estimés.

2. Procédé sur la revendication 1, caractérisé en ce que les paramètres a du modèle auto-régressif: P z(n) = Z aiy(n-i) = atyp(n) i=l avec at - ap, ap_.,...., a1] yp(n)t = [y(n-p),..., y(n-1)] o at est le vecteur des coefficients du modèle dont

l'ordre est p. o b(n) représente le processus généra-

teur du signal, o z(n) désigne les échantillons de la fonction représentant le phénomène, et o y(n) désigne les échantillons observes, sont régularisés par mise en oeuvre d'une contrainte consistant à minimiser un critère de la forme: 32(a) - Jj(a) + La) o J1(a) est la somme des moindres carrés [y(n)-z(n)]2, et o JIa) est une fonctionnelle régularisante tenant

compte de la connaissance a priori de la nature du phé-

nomène.

3. Procédé selon la revendication 2, caractérisé en ce que n(a) est de la forme: _(a) = " la - a012 o ao0 est une valeur prédéterminée, choisie a priori nulle ou constituée par le résultat obtenu pour la

tranche précédente.

4. Procédé selon la revendication 3, caractérisé

en ce qu'on prend " - p. es ou me est une valeur prédé-

terminée et "' la puissance estimée du signal.

,5. Procédé selon la revendication 3 ou 4, carac-

térisé en ce que la densité spectrale de puissance est calculée par la relation: p

g(f) - ô'/ I1l-Z a exp (-2infk)j 2].

6. Procédé selon la revendication 3, 4 ou 5,

caractérisé en ce qu'on traite simultanément les compo-

santes en phase et en quadrature d'un signal à bande

étroite démodulé.

7. Procédé selon la' revendication 3, 4, 5 ou 6, caractérisé en ce que, les signaux étant des signaux Doppler démodulés représentatifs de l'écoulement sanguin présentant un spectre de fréquence de quelques kHz, on utilise des tranches d'environ 64 échantillons prélevés

& une fréquence d'environ 40 kHz et on traite simultané-

ment les composantes en phase et en quadrature du signal démodulé.