text: Traitement numérique des signaux et identification de systèmes RMN ; Conception et développement d'un imageur pour l'IRM dédiée à bas champ Aktham Asfour
  To cite this version: Aktham Asfour. Traitement numérique des signaux et identification de systèmes RMN ; Conception et développement d'un imageur pour l'IRM dédiée à bas champ. Sciences de l'ingénieur [physics]. Université Joseph-Fourier - Grenoble I, 2002. Français. tel-00473212
  HAL Id: tel-00473212 https://tel.archives-ouvertes.fr/tel-00473212 Submitted on 14 Apr 2010
  
   UNIVERSITE JOSEPH FOURIER- GRENOBLE 1 N° attribué par la bibliothèque
  LLLLLLLLLLLI
  THESE pour obtenir le grade de
  DOCTEUR DE L'UNIVERSITE JOSEPH FOURIER Discipline: Sciences de l'ingénieur : Signal, Image préparée au laboratoire d'Electrotechnique de Grenoble dans le cadre de l'Ecole Doctorale
  « Electronique, Electrotechnique, Automatique, Télécommunications, Signal » présentée et soutenue publiquement
  par
  Aktham ASFOUR
  le 20 Juin 2002
  Titre:
  Traitement numérique des signaux et identification de systèmes RMN ; Conception et développement d'un imageur pour l'IRM dédiée à bas champ
  Directeurs de thèse: Kosai RAOOF et Jean Marc FOURN IER
  JURY
  M. Jean-Louis LACOUME M. André CONSTANTINESCO M. François GUILLET M. Marcel LOCATELLI M. Kosai RAOOF M. Jean-Marc FOURNIER
  Président Rapporteur Rapporteur Invité Directeur de thèse Directeur de thèse A la mémoire de mon père. A ma mère et à mes frères  Remerciements Le travail présenté dans ce manuscrit a été réalisé au sem du Laboratoire d'Electrotechnique de Grenoble (LEG) en collaboration avec le Laboratoire des Images et des Signaux (LIS).
  Je tiens à remercier chaleureusement Monsieur Kosai Raoof, Maître de conférence associé à l'Université Joseph Fourier de Grenoble, pour la confiance qu'il m'a témoignée en m'accueillant sous sa direction. Ses compétences scientifiques, sa disponibilité et ses qualités humaines furent primordiales au bon déroulement et à l'aboutissement de ce travail. Qu'il soit assuré de ma profonde gratitude et de ma sincère reconnaissance.
  Monsieur Jean-Marc Fournier, Professeur à l'Université Joseph Fourier de Grenoble a co-encadré cette thèse avec enthousiasme et bonne humeur. Je le remercie vivement pour la confiance qu'il m ' a accordée tout au long de cette thèse. Ses remarques pertinentes, les nombreuses discussions et ses qualités pédagogiques m'ont été très précieuses et encourageantes. Qu' il trouve ici l' assurance de mon profond respect et sympathie.
  J' adresse mes remerciement à : Monsieur Jean-Louis Lacoume, Professeur à l'Institut National Polytechnique de Grenoble, pour l'honneur qu'il m'a fait en assurant la présidence du jury de cette thèse.
  Monsieur André Constantinesco, Professeur à l'Université Louis Pasteur de Strasbourg, pour avoir accepté d'être rapporteur de cette thèse et d'avoir mis à profil sa longue expérience dans l'IRM pour l'évaluation de ce travail.
  Monsieur François Guillet, Maître de conférence à l'Université Jean Monnet de Saint Etienne, d'avoir accepté d'être rapporteur et pour avoir mobilisé son temps et ses compétences pour l'appréciation de ce travail.
  Monsieur Marcel Locatelli, Ingénieur au CEA de Grenoble, pour l'intérêt qu'il a manifesté à ce travail et pour avoir accepté d'être membre du jury. Je remercie chaleureusement : Les services électronique, mécanique et informatique du LEG pour leur aide précieuse, leur efficacité et leur sympathie. Une pensée particulière pour Bruno Mallet, qu'il soit assuré de ma profonde amitié. Le personnel administratif du LEG pour sa disponibilité, sa sympathie et sa bienveillance. Tous les membres du LEG pour l'ambiance chaleureuse qu'ils ont créée au laboratoire durant ces trois années. Tous les membres du LIS pour leur accueil lors de mes nombreux passages au LIS .
  Je suis très reconnaissant à Monsieur Jean-Louis Leviel et à tout le personnel du département Mesures Physiques de l'IUTl de Grenoble pour la confiance qu'ils m'ont accordée en m'accueillant en poste ATER durant la dernière année de cette thèse. Merci également au service du personnel de l'IUT 1 de Grenoble pour son aide et sa compréhension dans la réalisation des démarches administratives.
  Je remercie aussi Gilles Rostaing de m'avoir fait profiter de ses compétences et qualités pédagogiques en me confiant des tâches d'enseignements à ses côtés.
  Un grand merci à Makdam pour son aide précieuse lors de la rédaction de ce manuscrit et pour ses encouragements permanents.
  Je n'oublierai pas tous mes amis et compagnons au LEG que j'ai connus durant ces trois années et avec qui j'ai pu passer des moments agréables. Je leur souhaite bonne continuation et courage.
  Je ne saurai oublier tous mes amis qui se reconnaîtront pour leur présence, leur soutien et leur grande gentillesse. Je ne les citerai pas de peur d'en oublier quelqu'un.
  Je réserve mon dernier mot de remerciement pour ma mère et mes frères qui sont présents depuis toujours et à tout moment dans mon coeur.
  Merci à celles et à ceux que je n'ai pas cité et qui ont contribué de près ou de loin à ce travail. Table des Matières PRÉAMBULE 7 INTRODUCTION 9 LA RÉSONANCE MAGNÉTIQUE NUCLÉAIRE RMN : RÉALISATION D'UN SIMULATEUR DE SIGNAUX ••••••••••. 12 INTRODUCTION         . 13 DESCRIPTION DU PHÉNOMÈNE ET GÉNÉRALITÉS       13
  1. II.
  II. 1. II. 2.
  III.
  Notion de la fréquence de Larmor             13 Le phénomène de la résonance magnétique nucléaire (RMN)  15
  Bo
  II.2.1.
  Effet du champ magnétique statique
  II.2.2. II .2.3.
  Résonance magnétique: action d'un champ magnétique radiofréquence (RF)    16 La période de la précession libre et le signal RMN:      18
  sur un échantillon     15
  LES ÉQUATIONS DE BLOCH        20
  Il1.1 .
  Excitation radiofréquence RF                  20
  III.I .I. IILI.2.
  Caractéristiques de l'impulsion  20 Le repère tournant                      . 21
  Réponse du système à une impulsion RF.' l'équation de Bloch        22
  Il1.2. III.2.1 . IIL2.2. III.2J. III. 2.4.
  Comportement de l'aimantation macroscopique pendant l'impulsion d'excitation  22 Excitation On-resonance         23 Excitation Off-resonance                . 25 La période de la précession libre, l'aimantation mesurable  26
  IV. Du SIGNAL À L'IMAGE, LOCALISATION SPATIALE ET GRADIENTS DE CHAMP MAGNÉTIQUE   28 V. RÉALISATION D'UN SIMULATEUR DE SIGNAUX RMN    30
  V 1. V2. V3.
  Objectif et intérêt                          30 Principe et méthode de simulation     31 Résultats de simulation.' exemples de signaux RMN         34
  V.3.1 . V.3.2.
  Simulation d'une expérience avec une impulsion rectangulaire:  34 Simulation d'une expérience avec une impulsion de forme quelconque:  36
  DESCRIPTION ET CARACTÉRISTIQUES DE L'IMAGEUR DÉVELOPPÉ POUR L'IRM DÉDIÉE À BAS CHAMP 38 1. II.
  INTRODUCTION        39 DESCRIPTION DE L'IMAGEUR NUMÉRIQUE    39
  Il.1. II.2. II. 3.
  L'électroaimant                  . 41 L'antenne et le duplexeur   41 L'amplificateur RF de puissance  45
  IIJ.l. 11.3.2. II.3.3.
  Il.4. Il.5.
  La Carte DSP (Digital Signal Processor)   51 Le synthétiseur numérique      52
  11.5.1. II.5.2.
  Il.6.
  Le préamplificateur   46 L'étage sélectif       47 Caractéristiques mesurées de l'amplificateur total  48
  Intérêt du synthétiseur numérique  52 Le fonctionnement du synthétiseur AD9852  53
  Le récepteur numérique  56
  II.6.1. Rappel: récepteur analogique         56 11.6.2. Intégration d'un récepteur numérique     58 II .6.2.1. Description générale    58 n.6.2.2. Caractéristiques principales du récepteur numérique utilisé  59 11.6.2.3. Evaluation expérimentale des caractéristiques du récepteur     61
  RÉSULTATS: SIGNAUX RMN ET PREMIÈRES IMAGES 63 1. II.
  INTRODUCTION      64 SIGNAUX RMN ACQUIS AVEC L'IMAGEURNUMÉRlQUE  64
  II. 1. Il.2.
  Signal de décroissance d'induction FID  64 Les signaux d'écho   65
  II.2.1.
  Signal d'écho de spin    65 II.2.2. III.
  Signal d'écho de gradient.           68
  CONSTRUCTION DE L'IMAGE EN IRM    70
  IIU. Ill. 2.
  Sélection d'une coupe dans l'échantillon :  70 Codage de fréquence et codage par la phase: reconstruction de l'image  72
  II1.2.1 Codage par la fréquence        72 Codage par la phase     74 II1.2.2 m.2.3 Association du codage par la fréquence et du codage par la phase  74 III. 2.4. Acquisition de l'image, balayage du plan de Fourier:    75 llI.2.4.1. Principe du balayage             76 Ill. 2.4. 2. Considérations pratiques, paramètres d 'acquisition          77
  II!. 3. Les premières images obtenues avec le système numérique d'IRM par la séquence d'écho de gradient:      80 m.3.l. III.3.2.
  Réglage de l'angle de basculement de l'impulsion      80 La séquence d'écho de gradient              81
  IV. QUELQUES TRAVAUX EN COURS   83
  IDENTIFICATION DES SYSTÈMES RMN PAR FILTRAGE ADAPTATIF NON LINÉAIRE: MODÈLE DE SÉRIE DE VOLTERRA 85 1. II. III.
  INTRODUCTION  86 NON LINÉARITÉ DE SYSTÈMES RMN ET NÉCESSITÉ DE L'IDENTIFICATION NON LINÉAIRE    87 LA SÉRIE DE VOLTERRA APPLIQUÉE À L'IDENTIFICATION DE SYSTÈMES RMN  88
  IIU. III.2.
  Forme discrète de la série de Volterra:  88 Identification par série de Volterra de second ordre  89
  III.2.I. Méthode et schéma gériéral d'identification  89 III.2.2. Algorithme d'adaptation      91 m.2.2.1. Choix des pas d'adaptation pour les signaux RMN  93 111.2.2.2. Choix de la mémoire du système N  94 m.2.3 . Résultats de l'identification et discussions  95 I11.2.3 .1. Résultats  95 III.2.3 .2. Discussions   100
  Ill. 3. IV.
  Identification par série de Volterra de
  i
  me
  ordre    101
  MODÈLE DU FILTRE À RÉPONSE IMPULSIONNELLE INFINIE NON LINÉAIRE  105
  MODÉLISATION DES SYSTÈMES RMN PAR FILTRE À RÉPONSE lM PULSIONNELLE INFINIE (RII) NON LINÉAIRE: VALIDATION SUR DES SIGNAUX RÉELS 108 1 Il
  INTRODUCTION    109 MODÈLE DU FILTRE RH NON LINÉAIRE ET MÉTHODE D'IDENTIFICATION  . 109
  II.1 Il. 2
  Il. 3
  Description du modèle   109 Les approches du filtrage adaptatif RII non linéaire    109 L 'algorithme d'adaptation    111
  II.3.1 Les équations de mise à jour  II.3.2 Conditions de stabilité et pas d'adaptation variables  II.3.2. 1 Conditions de stabilité sur les pas d'adaptation  Il.3.2.2 Les pas d'adaptation variables pour les signaux RMN 
  III
  III 113 113 115 VALIDATION DU MODÈLE SUR DES SIGNAUX RÉELS RMN  116
  111.1 Ill. 2
  Exemple de signaux réels  116 Résultats d'identification  117
  CONCLUSIONS ET PERSPECTIVES  123
  ANNEXE  126
  BmLIOGRAPHIE  132 Préambule
  7 Le travail présenté ici a été réalisé au sein du Laboratoire d'Electrontechnique de Grenoble (LEG) et en étroite collaboration avec le Laboratoire des Images et des Signaux (LIS) pour la réalisation d'un imageur pour l'IRM à bas champ. Le projet a été partiellement supporté par le contrat européen QTEP ACK 1 pour le contrôle de qualité agroalimentaire.
  Le développement de l'imageur numérique n'aurait pas pu être accompli sans l'aide précieuse du Professeur André Constantinesco du Laboratoire de Biomécanique de CHU de Strasbourg qui n'a pas hésité à mettre à notre disposition la source de champ magnétique.
  1 QTEPACK
  QLKl- 2000-00936.
  8 Introduction
  9 Depuis sa découverte en 1946, la résonance magnétique nucléaire (RMN) n'a cessé de prouver ses qualités analytiques remarquables. Cette technique a permis à la spectroscopie de conquérir de nouveaux et de passionnants champs d'application. La RMN est devenue une des techniques analytiques les plus prospères pour un large domaine d'applications, s'étendant aujourd'hui de la physique de solide à toutes les branches de la chimie, la biologie moléculaire, le contrôle agroalimentaire non destructif et l'imagerie médicale. Le principe de cette technique est fondé sur la mesure de signaux en provenance de certains noyaux atomiques d'un échantillon en réponse à urie impulsion d'excitation électromagnétique radiofréquence. Dans le domaine médical, l'imagerie par résonance magnétique (IRM) produit des images du corps humain et des systèmes biologiques révélant la structure, le métabolisme, et la fonction des tissus et des organes internes. Son caractère non invasif et non ionisant associé à sa capacité de produire des images anatomiques avec une qualité sans précédente l'ont rendu plus attirante et plus fascinante que beaucoup d' autres techniques d'imagerie. Cependant, cette technique reste aujourd'hui réservée aux grands départements de radiologie et ses indications sont davantage restreintes à certaines pathologies. Cette limitation de disponibilité est surtout liée au coût très élevé des équipements dû aux aimants, souvent supraconducteurs produisant les moyens et les hauts champs magnétiques utilisés, ainsi qu'à l'électronique d'acquisition, qui est jusqu'à présent, en partie analogique. Ces imageurs «corps entiers » ne sont pas, en plus, bien adaptés à l'imagerie des petits organes (poignet, main, genou  ). L'IRM dédiée à bas champ constitue une alternative peu onéreuse et permet donc d'assurer une meilleure disponibilité de l'IRM dans ces cas. Ce coût peu élevé peut également rendre cette technique très attractive et accessible dans des domaines autres que médicaux tel que le contrôle agroalimentaire non destructif. L'utilisation de l'IRM à bas champ est, bien évidemment, synonyme d'un faible rapport signal sur bruit (RSB), et par conséquent d'une qualité d'image inférieure pour être exploitée dans les applications de diagnostic ou de contrôle. Le travail présenté dans ce manuscrit porte, dans sa première partie, sur le développement d'un système numérique d'IRM à bas champ permettant l'amélioration du rapport signal sur bruit et de la stabilité du système. Le choix d'une solution numérique est argumenté par le besoin d'élaborer un imageur peu coûteux et de faible encombrement pouvant être facilement intégré dans tout environnement aussi bien médical, tel qu'un cabinet de médecin, qu'industriel comme le contrôle agroalimentaire.
  10 L'originalité du système développé réside principalement dans l'utilisation d'un synthétiseur numérique de fréquence, avec contrôle absolu de phase, et d'un récepteur numérique qui constituent les deux éléments clés d'un imageur. Ceux-ci mettent à profit le très bon rapport qualité-prix des systèmes numériques haute fréquence autorisé par les développements technologiques incessants des télécommunications et favorisé par la demande croissante dans ce domaine.
  Parallèlement à cette phase de conception et de développement instrumental, nous avons conduit, en se basant sur les techniques adaptatives en traitement du signal, une étude d'identification des systèmes RMN. L'objectif étant de permettre la prise en compte de la non linéarité de tels systèmes dans le traitement et l'analyse, jusqu'à présent linéaires, du signal RMN. Le travail reporté dans ce document peut alors être organisé en deux parties. La première partie est consacrée au développement de l'imageur numérique. Dans cette partie, seront rappelés, dans le premier chapitre, les principes physiques fondamentaux du phénomène de la résonance magnétique nucléaire qui servent de support pour toute personne abordant l'IRM. Ces fondements nous ont été nécessaires pour la mise au point et l'implantation d'un simulateur de signaux RMN. Dans le deuxième chapitre, nous exposons la conception de l'imageur et nous détaillons les solutions instrumentales numériques mises en oeuvre et l'évaluation de leurs caractéristiques de performances. Des exemples de signaux RMN acquis, la technique de reconstruction de l'image ainsi que les premières images réalisées avec cet imageur sont donnés dans le troisième chapitre. La deuxième partie, constituée des chapitres IV et V, traite la problématique de l'identification non linéaire des systèmes RMN à travers des techniques adaptatives de traitement du signal. Ainsi, nous introduisons, dans le quatrième chapitre, un modèle non linéaire basé sur la représentation d'un système RMN par un filtre de Volterra et nous le comparons au modèle d'un filtre non linéaire à réponse impulsionnelle infinie (RH). A l'issu cette comparaison, le filtre RH est retenu. Sa description détaillée et sa validation sur des signaux RMN réels font l'objet du cinquième chapitre.
  11 Chapitre 1
  La résonance magnétique nucléaire RMN : Réalisation d'un simulateur de signaux
  12 1.
  Introduction Le travail présenté tout au long de ce manuscrit consiste dans le développement
  instrumental et le traitement du signal pour l'IRM. Il est alors nécessaire d'introduire les principaux fondements théoriques indispensables à la compréhension du mécanisme de génération des signaux RMN. Ces fondements serviront, en particulier, dans ce chapitre à la mise en place d'un simulateur de signaux. Nous allons donc, dans un premier temps, et après une description succincte du phénomène physique de la RMN, présenter les équations régissant ce phénomène, c'est-à-dire les équations de Bloch, et rappeler leurs solutions dans des cas particuliers. Dans un deuxième temps, nous décrivons le simulateur de signaux dont les principes sont basés sur la résolution numérique des équations de Bloch. La validité du simulateur sera exposée à travers quelques exemples de signaux simulés.
  Il.
  Description du phénomène et généralités Il.1.
  Notion de la fréquence de Larmor
  Il est bien connu que plusieurs noyaux d'intérêt biologique possèdent un moment magnétique. Il s'agit des noyaux d'hydrogène (IH) formé d'un proton, de l'isotope de carbone l3C, d'azote
  eN), de fluor 4
  ( 19p ),
  de sodium
  e 3
  N) et de phosphore
  elp).
  Selon la mécanique classique, le moment magnétique de tels noyaux, lorsqu'ils sont placés dans un champ magnétique statique
  Bo,
  présente, dans un référentiel fixe de
  laboratoire (OXYZ), un mouvement de rotation appelé précession. Ceci est illustré dans figure 1.1 où le moment magnétique, désigné
  ïJ, balaye un cône d'axe parallèle au vecteur champ
  magnétique Bo, et ce avec une fréquence caractéristique appelée fréquence de Larmor [Guinet, 1992].
  13 y K
  J-1 Figure 1.1
  Précession d'un moment magnétique j1 d'un noyau placé dans un champ
  magnétique statique Bo.
  Cette fréquence est donnée par la relation fondamentale dite de Larmor:
  (1.1 )
  où y est le rapport gyromagnétique du noyau et Bo est le module du vecteur champ magnétique. Cette relation est parfois donnée en fonction de la vitesse angulaire de précession
  (1.2)
  où le signe (-) indique que la précession du moment s'effectue dans le sens de rotation des aiguilles d'une montre. En examinant la figure 1.1, nous constatons que, lors de la précession, la proj ection Jlz du moment magnétique sur l'axe Z reste constante, alors que la projection transversale Jl xy dans le plan XOY tourne dans le sens inverse au sens trigonométrique à la vitesse
  angulaire
  (00.
  14 Il est alors possible de caractériser un noyau par la fréquence de Lannor
  Jo
  de son
  moment magnétique, en prenant soin de préciser la valeur du champ magnétique statique Bo' On se limitera par la suite, et sans perte de généralité, au cas du noyau d'hydrogène (un seul proton) dont le rapport gyromagnétique
  r vaut 267.53 x 10 6 rad. S*I . TI. Pour une valeur
  de Bo égale à 0.1 T (Tesla), la fréquence de Larmor du proton vaut donc 4.26MHz.
  Il.2.
  Le phénomène de la résonance magnétique nucléaire (RMN)
  Il.2.1. Effet du champ magnétique statique
  En l'absence du champ magnétique externe
  Bo
  sur un échantillon
  Bo, les moments magnétiques des noyaux
  au sein d'un échantillon sont orientés de manière aléatoire dans l'espace si bien que leur somme est nulle: aucune aimantation macroscopique n'est, par conséquent, observable au voisinage de l'échantillon. Placée dans un champ magnétique statique homogène, la population de moments magnétiques des noyaux de cet échantillon sera le siège de l'apparition d'une aimantation macroscopique
  M
  parallèle et de même sens que
  longitudinale celle donnée par la direction de
  Bo,
  Bo.
  En prenant pour direction dite
  on dira que
  M
  n'a pas de composante
  transversale, i.e. sa projection dans le planXOY est nulle [Liang, 1999]. L'apparition de cette aimantation macroscopique
  M est attribuée au fait que dans le
  mouvement de précession, les moments magnétiques des protons sont répartis en deux populations : l'une où les moments
  ï1
  ont leurs projections longitudinales f.iz parallèles à
  Bo.
  C'est le cas de la figure LI. On dit que ces moments ont un niveau de basse énergie. l'autre où les moments
  ï1
  ont leurs projections longitudinales f.iz antiparallèles à
  Bo. Cette catégorie de moments se trouve dans un niveau de haute énergie. La répartition des moments dans le sens parallèle et antiparallèle n'est pas égale : les moments magnétiques parallèles à
  Bo
  sont légèrement majoritaires. Bien qu'infime, cette
  15 différence est à l'origine du vecteur d'aimantation macroscopique M de même sens que Bo [Nguyen, 1996] L'absence d'une composante transversale M xy de cette aimantation est due au fait que les composantes f-l XY des moments magnétiques j1 sont déphasées aléatoirement dans le plan XOY. Ainsi, l'aimantation macroscopique induite par
  13
  0
  se réduit à une composante
  longitudinale M z due aux composantes f-l z des moments magnétiques des protons en surnombre dans le niveau de basse énergie (état parallèle). Dans cette situation, on dit que l'échantillon est à l'équilibre, et nous désignons par M~ la valeur de cette composante qui est proportionnelle à la densité, p , de protons dans l'échantillon. L'amplitude de
  M
  étant trop faible pour être mesurée à l'équilibre, des méthodes
  d' excitation perturbant l'équilibre du système sont alors nécessaires pour détecter cette aimantation. On parlera donc de résonance magnétique.
  Il.2.2. Résonance magnétique: action d'un champ magnétique radiofréquence (RF) Le principe de la résonance magnétique protonique est d'agir spécifiquement sur les moments magnétiques des noyaux d'hydrogène à l'aide d'impulsions de champ magnétique tournant BI [Liang, 1999] [Guinet 1992]. Ce champ doit être dans un plan normal au champ statique
  130 , Son amplitude est de plusieurs ordres de grandeur plus faible que celle de 130
  tourne autour de ce dernier avec une vitesse angulaire
  OJ RF •
  et
  Il s'agit d'un champ appliqué de
  manière transitoire et appartenant à la gamme des radiofréquences, d'où le nom d'impulsion radiofréquence (RF) régulièrement employé pour désigner BI . L'action d'une impulsion RF de fréquence égale à ou proche de la fréquence de Larmor, se traduit par la décroissance de l'aimantation longitudinale M z et l'apparition d'une composante transversale de l'aimantation M XY tournant autour de OJ o
  comme le montre la figure I.2. Le mouvement de l'aimantation
  est une précession-nutation: elle tourne autour de entre
  M et 130
  130
  à la fréquence
  130
  à la vitesse angulaire
  M , pendant l'excitation, OJ o '
  tandis que l'angle a
  augmente proportionrtellement à la durée de l'impulsion" suivant la relation
  a = yB I " • L'extrémité du vecteur
  M décrit ainsi une spirale sur une surface sphérique. 16 z
  y
  Figure 1.2
  Evolution de l'aimantation macroscopique durant l'application de l'impulsion d'excitation.
  L'angle a, caractérisant l'impulsion, est appelé angle de basculement. Les impulsions sont souvent désignées par leur angle de basculement. Ainsi, on parle d'impulsion 90° si BI et
  'f
  sont choisis de telle sorte que a soit égal à 90°. A la fin d'une telle impulsion, on trouve
  uniquement une composante transversale d'aimantation. L'impulsion 180° correspond au retournement du vecteur
  111 qui, en fin d'impulsion, devient parallèle et opposé
  à Bo , soit
  donc une aimantation transversale nulle [Guinet, 1992]. L'évolution de l'aimantation longitudinale M z sous l'action du champ magnétique tournant peut s'expliquer par le fait que l'impulsion RF fournit, durant son application, une énergie nécessaire pour permettre la transition des moments magnétiques des noyaux de l'état parallèle (basse énergie) vers l'état antiparallèle (haute énergie) [Liang, 1999]. Lorsque la moitié des moments en surnombre est passé à l'état antiparallèle, les deux populations parallèle et antiparallèle sont en égalité et la composante M z s'annule: une impulsion 90° a été donc appliquée. Lorsque tous les moments en surnombre sont orientés antiparallèlement, il y a inversion de la composante longitudinale, et la composante transversale est donc nulle : on dit qu'une impulsion 180° a été appliquée.
  17 L'évolution de l'aimantation longitudinale ainsi présentée ne permet pas de décrire celle de la composante transversale M xy . Pour expliquer la naissance de cette dernière, on doit se rappeler qu'à l'équilibre, les composantes f.1 xy des moments des noyaux en précession dans le plan transversal sont déphasées aléatoirement. La résonance magnétique peut être vue comme le mécanisme d'établissement d'une cohérence de phase au sien de ce groupe de moments magnétiques [Liang, 1999]. En effet, lorsque les moments changent d'état d'énergie, ils le font en se mettant en phase comme si leurs projections transversales étaient poussées à l'unisson à leur fréquence de précession par une force extérieure. Cette force est apportée, bien évidemment, par le champ RF tournant de la même manière que *les composantes transversales qui se trouveront aussitôt en phase donnant ainsi naissance à la composante transversale M xy .
  Il.2.3. La période de la précession libre et le signal RMN:
  Après un basculement de l'aimantation
  M
  par une impulsion RF, l'échantillon est
  placé dans un état hors équilibre. Le processus de retour à l'équilibre se caractérise par une précession libr/ de l'aimantation
  M
  autour du champ
  Bo
  à la fréquence co o; une
  augmentation de l'aimantation longitudinale M z pour tendre vers sa valeur d'équilibre M~, appelée relaxation longitudinale,' et une destruction de l'aimantation transversale M xy , appelée relaxation transversale. Les processus de relaxation sont quantifiés par deux constantes de temps
  ~
  et T2 dites de relaxation longitudinale et de relaxation transversale,
  respectivement. Ces sont les temps au bout desquels M z et M xy ont atteint 63% de leurs valeurs à l'équilibre. A titre d'exemple, pour des tissus biologiques,
  ~
  peut varier de
  quelques 100ms à quelques secondes, alors que T2 est environ dix fois plus petit que
  ~
  [Liang, 1999] [Nguyen, 1996]. Lors de la relaxation, la composante transversale M xy décrit, dans l'" ulan XOY, un mouvement de rotation à la fréquence de précession
  Jo
  accompagné par la décroissance de
  son amplitude. Ce mouvement provoque, grâce au phénomène d'induction magnétique, l'apparition dans une antenne de réception, placée dans le plan XOY, d'un signal appelé signal de décroissance d'induction libre ou FID (Free Induction Decay en anglais) (figure 1.3). Ce
  1 Précession
  libre désigne la précession de l'aimantation
  M en l'absence de l'impulsion d'excitation RF.
  18 signal qui s'apparente à une sinusoïde amortie et dont la fréquence est celle de la précession est généralement détecté en quadrature (figure 1.3). Le signal résultant est l'image du phénomène de relaxation et décrit une fonction exponentielle amortie d'une constante de temps T2 et dont l'amplitude maximale est proportionnelle à M~, soit donc à la densité protonique p .
  , Z
  M
  y
  x Antenne réceptrice
  1
  Détection en quadrature
  Décroissance en
  e
  2
  T
  c=====> t
  Signal de décroissance d'induction (FID) à la fréquence
  Enveloppe du signal
  fa
  Figure 1.3 Trajectoire de l'aimantation M dans le repère fixe pendant la période de la précession libre et le signal FID détecté.
  En pratique, la décroissance exponentielle du signal s'effectue avec une constante de temps effective T2• plus petite que T2 et ce en raison des inhomogénéités du champ
  19 magnétique statique
  Ba
  qui accélèrent la décroissance de la composante transversale [Guinet,
  1992].
  III.
  Les équations de Bloch Dans la description de la physique, l'évolution de l'aimantation durant l'application de
  l'impulsion RF est régie par les équations de Bloch. C'est ce que nous détaillons dans ce paragraphe.
  111.1. Excitation radiofréquence RF 111.1.1. Caractéristiques de l'impulsion
  Nous venons de voir dans le paragraphe IL2.2 qu'une impulsion RF désigne l'application transitoire, de quelques microsecondes à quelques millisecondes, d'un champ magnétique
  BI (t)
  orthogonal à
  Ba'
  Ce champ, appelé encore excitation, est créé
  généralement par un courant sinusoïdal parcourant un fil conducteur. Dans le cas d'une bobine (solénoïde) d'axe parallèle à l'axe Ï du repère fixe ( voir figure LI), l'expression de
  BI
  est de la forme [Liang, 1999] :
  (1.3)
  où: Bt(t) est la fonction d'enveloppe de l'impulsion dont la dépendance par rapport au
  temps signifie que l'amplitude de ce champ est variable pendant sa durée d'application. On dit que l'impulsion est modulée en amplitude OJ RF
  est la fréquence de l'oscillation (fréquence de la porteuse)
  cp est la phase initiale
  Le champ BI donné par l'équation (L3) est dit polarisé rectilignement car il oscille linéairement le long de l'axe Ï. Ce champ peut se décomposer en deux composantes polarisées circulairement et tournantes dans deux sens opposés. Seule la composante tournante dans le sens des aiguilles d'une montre (sens de la précession de Larmor) est
  20 efficace pour agir sur le moment magnétique du noyau, et ce à condition que la fréquence soit égale où proche de la fréquence de Larmor [Liang 2000]. Cette composante
  OJ RF
  s'exprime par la relation2 [Ernst, 1989] [Liang 1999] :
  (1.4)
  D'une façon générale, les d'impulsions sont désignées, en plus de leur angle de basculement, par leur fonction d'enveloppe
  Ble (t).
  Parmi les fonctions les plus répandues,
  nous citons la fonction rectangulaire (porte) et la fonction sinus cardinal (sine). Nous parlerons donc, à titre d'exemple, d'impulsion rectangulaire 180° ou d'impulsion sine 90°.
  111.1.2. Le repère tournant Il est commode d'introduire, afin de comprendre plus facilement l'action du champ magnétique RF sur l'aimantation, un repère oxyz appelé repère tournant (figure lA). C'est un système de coordonnées dont le plan transversal xoy tourne dans le sens inverse du sens trigonométrique à la vitesse angulaire
  OJ RF
  [Ernst, 1989]. L'axe z de ce repère est confondu
  avec l'axe Z du repère fixe alors que les axes x et y tournent donc dans le plan XOY comme
  BI'
  Ainsi, dans ce repère tournant, le vecteur
  BI
  a une direction fixe. Ses deux composantes
  s'écrivent:
  2
  BI (t)
  Bt (t)
  = Bt (t) cos(cp)
  (1.5)
  Bi(t)
  = Bt(t)sin(cp)
  (1.6)
  désigne désormais la composante utile du champ magnétique radiofréquence.
  21 y
  Repère de laboratoire OXYZ et repère tournant oxyz.
  Figure 1.4
  La phase cp de l'impulsion RF est définie par l'excursion du vecteur
  BI
  de l'axe x par
  rapport à l'axe y. Ce paramètre est d'une importance cruciale dans les expériences RMN et doit être contrôlé adéquatement. C'est ce que nous verrons dans le chapitre II.
  111.2. Réponse du système à une impulsion RF : l'équation de Bloch 111.2.1. Comportement de l'impulsion d'excitation
  l'aimantation
  macroscopique
  pendant
  Dans le repère tournant, l'équation de Bloch décrivant l'évolution de l'aimantation macroscopique durant l'application de l'impulsion se présente sous la forme vectorielle suivante [Liang, 1999] :
  dM /ollr _ _ _ =y(M/ollr
  ~
  x
  (M ~ollr i + My/our ]) BejJ ) + ----::.-~
  22
  /M /our
  \'
  z
  -
  ~
  M
  o z
  )k (1.7) où
  M
  lou
  ,
  = (MIOJ/' My x '
  lou
  ,
  ,
  MIOJ/,) est le vecteur d'aimantation dans le repère tournant et Ëeff est z
  le champ magnétique effectif. C'est, en effet, le champ total vu par l'aimantation dans le repère tournant. Son expression est la suivante [Liang, 1999]:
  (1.8)
  Le phénomène de relaxation peut être négligé pendë.Ilt l'application de l'impulsion. La durée de celle-ci étant généralement très petite par rapport à
  T
  et à T2 et l'équation de Bloch
  se réduit à :
  dM
  --
  -dt= y(MxB ,tr )
  (1.9)
  e))
  où
  M
  désigne maintenant le vecteur aimantation dans le repère tournant. Dans cette dernière
  équation, nous avons fait abstraction de l'indice tour car nous travaillons désormais dans le repère tournant. L' équation de Bloch n'a en général de solution analytique que dans quelques cas où la fonction d'enveloppe
  Bt (t)
  a une forme particulière. C'est le cas par exemple des impulsions
  rectangulaires. Dans le cas général, il conviendra donc de faire appel à des méthodes de résolution numérique. Dans ce qui suit, nous nous intéressons au comportement de l'aimantation sous l'action du champ RF dans le cas d'une impulsion rectangulaire. Pour cela, nous allons procéder à la résolution de l'équation de Bloch en distinguant deux cas selon que la fréquence de l'excitation est égale (On-resonance) ou légèrement différente (Off-resonance) de la fréquence de précession. 111.2.2. Excitation On-resonance
  Par souci de simplification, nous supposons que le champ BI est appliqué selon l'axe
  T du repère tournant ( cp = 0) pendant une durée
  r . Le terme On-resonance est employé pour
  indiquer que la fréquence de l'impulsion est parfaitement égale à la fréquence de précession
  23 ({1)RF
  = (1)0). Dans ce cas, le champ effectif (équation 1.8) se réduit au champ radiofréquence
  et s'écrit comme:
  (1.10)
  En prenant en compte les conditions initiales de l'équilibre, soit M x(0) = M /0) = 0 et Mz (O) = M~, il est possible de montrer que l'évolution temporelle des trois composantes de l'aimantation pendant l'excitation peut s'écrire sous la forme [Liang, 1999] :
  t
  M/t) = M~sin( fyBt(O)dO )
  (1.11)
  o t
  Mz (t) = M~cos( f yBt(O)dO ) o
  Dans le cas d' une impulsion rectangulaire d'amplitude Bt, ces équations deviennent: Mx (t) = 0
  = M~ sin(yBtt) Mit) = M~ cos(yBtt)
  M /tJ
  (1.12)
  Ces équation laissent apparaître que le mouvement de l'aimantation dans le repère tournant est une rotation (figure 1.5). A un instant t donné pendant l'impulsion, tout se passe comme si l'aimantation M~ avait tourné dans un plan perpendiculaire à BI' soit le plan yoz, d'un angle yBtt. A la fin de l'impulsion, l'aimantation aura alors tourné de l'angle de basculement a = yBt-r. Pour une impulsion de forme quelconque et de durée -r , cet angle est donné par :
  r
  (Ll3)
  a = y fBt(O)dO o
  24 z
  ~M(t) 
    
  a=yBtr
  yBtt.// .'  
  M(r)
  ./.'
  y
  x Figure 1.5 Evolution de l'aimantation dans le repère tournant pendant une excitation rectangulaire de durée r. L'angle de basculement a est défini comme le plus petit angle entre
  l'aimantation M(r) à la fin de l'impulsion et l'axe z.
  L'aimantation qui n'avait à l'équilibre qu'une seule composante longitudinale Mz
  = M~,
  possède à la fin de l'impulsion, en plus d'une composante longitudinale, une
  composante transversale M xy . Dirigée suivant l'axe y (My (r) = M xy (r) ), la composante M xy peut avoir en général une direction quelconque dans le plan xoy. Cette direction, qui représente en réalité la phase du signal RMN reçu, dépend bien évidemment de la direction du champ radiofréquence dans le repère tournant, soit donc de la phase de l'impulsion d'excitation.
  111.2.3. Excitation Off-resonance Dans les expériences de RMN, les excitations sont supposées de type On-resonance. Néanmoins, les inhomogénéités du champ statique rendent très difficile la réalisation de telles excitations dans tout le volume de l'échantillon à étudier. Les impulsions sont donc souvent de type Off-resonance, i.e.
  OJ RF :1;
  co o ' Soit ~OJo
  = CO o -
  OJ RF
  la quantité traduisant cet écart de
  fréquence, le champ effectif ËeJJ peut alors s'écrire sous la forme:
  25 -
  BefJ
  e -: o= -ÔCû - k + BI (t)l
  (1.14)
  y
  Comme nous l'avons dit, une solution analytique de l'équation de Bloch ne peut être obtenue dans le cas général. Cependant, on peut intuitivement prédire, par analogie avec l'excitation On-resonnance, qu'à la fin de l'impulsion, tout se passe comme si l'aimantation avait tourné dans un plan perpendiculaire à la direction du champ effectif. Cette direction étant celle de BI dans le cas de l'excitation on-resonance. Ainsi, la composante transversale de l'aimantation n'est pas forcément dirigée selon l'axe oy. Le terme
  ô010 = 01 0
  -
  Cû RF
  dans l'expression de BefJ est souvent négligé lorsque le
  champ radiofréquence est suffisamment intense (yBt
  » 1010 -01 1) RF
  [Ernst, 1989] [Canet,
  1991]. Avec cette approximation, l'étude se ramène au cas On-resonance traité ci-dessus.
  Pour plus de détails sur les équations de Bloch et leurs solutions, le lecteur intéressé pourra se référer à la littérature RMN [Liang, 1999] [Ernst, 1989].
  111.2.4. La période de la précession libre, l'aimantation mesurable
  Après la fin de l'impulsion d'excitation, le retour à l'état d'équilibre du système peut être également décrit par l'équation de Bloch. En fait, en l'absence du champ radiofréquence BI ' l'expression du champ effectif B efJ devient: (1.15)
  L'équation de Bloch conduit ainsi au système d'équations suivant [Liang, 1999]:
  (1.16)
  26 Si on considère maintenant que l'origine des temps, t=O, est pris à l'instant où l'excitation, supposée rectangulaire, s'annule et que le terme Off-resonnace est négligeable par rapport au champ radiofréquence, i.e. rBI
  » 1mo - m l, nous obtenons les conditions RF
  initiales suivantes :
  M x(0) = 0 MiO) = M~ sin(a)
  (1.17)
  MlO) = M~cos(a)
  On montre alors que la solution du système d ' équations (1.16) peut être donnée comme [Ernst, 1989] : -(
  - (
  Mit) = M ~ sin (a )sin( L\mot) e T2
  (1. 18) -(
  - (
  Mlt) =M~cos(a )e 1i +M~(1 -e 1i )
  En introduisant l'aimantation transversale complexe, soit M xy
  = Mx + iM y'
  nous
  obtenons: -(
  M xy (t)
  = M~ sin(a)e-
  i6 (j)o(
  e
  T2
  (I.19)
  Une expression analytique de l'aimantation mesurable est donc obtenue car nous nous sommes basés sur une impulsion rectangulaire. Dans le cas général, ceci nécessitera la mise en oeuvre de méthodes numériques de résolution. Le simulateur de signaux RMN que nous présentons plus loin de ce chapitre en constitue une. Il implante, en fait, des opérateurs matriciels pour calculer l'effet de l'impulsion et décrire l'évolution de l'aimantation pendant la période de précession libre.
  27 IV.
  Du signal à l'image, localisation spatiale et gradients de champ magnétique Nous venons de voir dans les paragraphes précédents, qu'une expérience RMN
  détennine la densité de protons et les temps de relaxation T2 et
  ~
  de l'ensemble de
  l'échantillon étudié. La formation d'une image (image de la densité de proton dans l'échantillon
  3
  )
  à partir
  d'un signal nécessite la localisation spatiale des signaux émis par tous les volumes élémentaires, ou voxels, de l'échantillon. Il s'agit donc de codage spatial [Ravier!, 1994]. La méthode la plus couramment utilisée pour le codage spatial est dite de transformée de Fourier. Elle est basée sur la relation de Larmor liant la fréquence de résonance, f , d' un voxel à l'intensité du champ magnétique statique B(X,Y ,Z) appliqué.
  f( X ,Y, Z)
  = yB(X,Y, Z)
  (1.20)
  où X, y, Z sont les coordonnées du vox el considéré dans le repère du laboratoire. Cette relation montre que, si B est unifonne et vaut Bo, tous les protons de l' échantillon résonnent à la même fréquence, f o' et ce quelle que soit leur position dans l' espace. Le principe de l'imagerie est d'imposer une variation spatiale sur B créée par des gradients de champ magnétique. Cette variation génère ainsi une variation spatiale de la fréquence de résonance des protons. Lorsque la variation spatiale de B est linéaire dans les trois directions de l'espace, la fréquence du signal émis par un vox el dépend linéairement de la position de ce vox el. D'une manière générale, pour produire un gradient de champ magnétique, on superpose au vecteur champ magnétique
  Bo
  un vecteur champ magnétique
  b parallèle à Bo.
  L'amplitude b du vecteur b , qui est de l'ordre du millitesla (mT), est beaucoup plus faible que Bo' Contrairement à Bo, b varie selon les trois axes X, Y,Z dans l'espace qui contient l'échantillon à imager
  On peut également obtenir des images dites ~ et T2 représentant respectivement l'image de la distribution du temps de relaxation longitudinale et du temps de relaxation transversale au sein de l'échantillon.
  3
  28 Pour fixer les idées, considérons l'exemple où b varie linéairement selon l'axe Y (figure 1.6) selon l'équation:
  (1.21)
  Pour tous les points d'un plan de cote Y, i.e. parallèle au plan XOZ, le vecteur champ magnétique a la même valeur B donnée par l'équation :
  (1.22)
  La constante G y est appelée gradient du champ magnétique selon Y. Elle vaut donc :
  G
  = aB
  (1.23)
  ay
  y
  Le gradient Gy peut être positif ou négatif et s'exprime habituellement en millites1as par mètre (mT/m) et Z
    . . . 
  y
  o
  ********rJlll****.J-----I--LLL~
  Figure 1.6
  Gradient de champ magnétique selon l'axe Y.
  29 En pratique, les gradients sont créés par des bobines appelées bobines de gradient. Elles sont généralement conçues pour avoir une valeur de b nulle à l'origine [Guinet, 1992]. Chaque instrument d'IRM possède trois gradients dans les trois directions de l'espace. Nous allons voir dans le troisième chapitre de ce manuscrit comment l'utilisation adéquate de tels gradients de champ permet de limiter la région excitée et de coder spatialement le signal reçu. On parle ainsi de séquence d'imagerie. En effet, une séquence d'imagerie est une suite d'impulsions RF et de gradients appliqués à un corps pour créer un nombre suffisant de signaux RMN nécessaire à la reconstruction d'une image. Généralement, une séquence peut être divisée en plusieurs cycles répétés avec un temps de répétition approprié TR • Le cycle commence, dans la plupart des cas, par une impulsion RF. Le nombre de répétitions du cycle dépend du type de séquence utilisée.
  v.
  Réalisation d'un simulateur de signaux RMN
  V.1.
  Objectif et intérêt
  La simulation en IRM est un complément important à l'expérimentation. Outre sa précieuse contribution dans la compréhension générale de la théorie de la résonance magnétique nucléaire, la simulation permet en réalité d'étudier et d'analyser séparément les paramètres d'influence pouvant intervenir dans l'expérience. Ces paramètres n'étant pas tous accessibles à l'utilisateur d'un appareil d'IRM [Brenner, 1997], [Summer, 1986]. La simulation trouve notamment tout son intérêt pour répondre aux diverses questions pouvant se poser pendant l'expérience [Bittoun, 1984]. Une liste non exhaustive de ces questions comprend les interrogations suivantes: Qu'est ce qu'on mesure par chaque séquence; quel type de distorsions d'image est provoqué par une modification des caractéristiques de l'impulsion RF ou par un changement dans la séquence des gradients ; lesquelles de ces distorsions sont causées par des contraintes technologiques et lesquelles sont dues à des raisons théoriques? La réponse à de telles questions nécessite la connaissance du comportement de l'aimantation régie par l'équation de Bloch dont la solution ne peut être établie que pour un nombre très limité de types d'excitation. La simulation permettra ainsi de s'affranchir de cette limitation.
  30 De plus, la spécification, l'optimisation des séquences d'imagerie voire encore le dimensionnement de nouveaux systèmes peuvent être établis à l'aide de la simulation qui pourra, d'ailleurs, constituer un outil pédagogique dans des environnements médicaux et techniques.
  Nous allons dans ce qui suit décrire la mise en place et la validation d'un simulateur de signaux RMN et nous présenterons et comparerons des exemples de signaux simulés avec des résultats reportés dans la littérature. Ces signaux seront utilisés pour la validation des approches d'identification non linéaire de systèmes RMN traitées dans le chapitre IV.
  V.2.
  Pri ncipe et méthode de simulation
  La simulation en RMN consiste à résoudre numériquement les équations de Bloch, dans le repère tournant, pour chaque point (vox el) de l'échantillon étudié et ce à chaque instant de la séquence. L'accumulation des signaux élémentaires correspondant à tous les voxels de l'échantillon simulé permet de calculer un signal RMN global [Bittoun, 1984]. Il nous faut donc d'abord définir l'échantillon en question. Il s'agit de définir un objet ID, 2D ou 3D pouvant, bien entendu, avoir la forme d'une sphère, d'un cylindre ou d'autres formes géométriques définies par une expression mathématique. Cet objet est discrétisé en petits éléments de volume (voxels) auxquels sont attribuées des caractéristiques essentielles. Celles-ci comprennent notamment les dimensions, le vecteur aimantation locale
  m, la valeur
  de l'aimantation à l'équilibre ma (proportionnelle à la densité de protons dans le voxel), et les temps de relaxation locaux 1; et T2 • La deuxième étape de la simulation consiste à résoudre l'équation de Bloch pour chaque voxel de l'objet ainsi défini et ce dans le cas d'une excitation quelconque. Pour cela, nous considérons le fait que la durée de l'impulsion RF est petite par rapport aux temps de relaxation longitudinale et transversale. Avec cette hypothèse, l'équation de Bloch relative à un vox el s'écrit, dans le repère tournant:
  dm dt où
  = y(mx Beff}
  (1.24)
  m désigne maintenant le vecteur aimantation du voxel.
  31 La méthode de résolution de cette équation fait appel à un opérateur matriciel décrivant la rotation, dans le repère tournant, de l' aimantation
  m sous l'effet de l'excitation.
  Soit m( X, Y, Z, t) le vecteur aimantation au point (X, y, Z) à un instant donné t. A la fin d'une impulsion rectangulaire de durée r, celui-ci s'exprime à l'aide de cet opérateur par la relation suivante [Bittoun 1984] [Brent, 1986], [Brenner 1997] :
  m(X,Y,Z ,t + r) = RRF(X,Y,Z), m(X,Y,Z,t)
  (1.25)
  L 'opérateur RRF ' appelé opérateur d'impulsion, est donné par la matrice suivante :
  COS
  2
  RRF(X, Y,Z) = [
  (j3) + sin 2 (fJ) cos(0)
  sin( 0) sin(j3)
  - sin(j3) cos(fJ)[I - cos(
  0)]]
  - sin(o)sin(j3)
  cos(o)
  - sin(0) cos(fJ)
  - sin(j3) cos(j3) [1 - cos( 0)]
  sin(0) cos(j3)
  sin 2 (j3) + cos 2 (fJ) cos(0) (1.26)
  où 0 et j3 sont des angles qui dépendent du champ effectif
  Bef!
  et que nous explicitons dans
  ce qui suit. Nous reprenons l'expression du champ effectif donnée par l'équation 1.8 tout en prenant en compte que, lors d'une séquence d'imagerie, le champ magnétique statique en un point (X, y, Z) a la forme générale suivante:
  où Mo caractérise l'inhomogénéité du champ statique principal
  Bo
  en ce point; Gx ' Gy et
  Gz sont les trois gradients appliqués pendant la séquence d'imagerie pour coder l'espace. Le champ effectif
  Bef!
  en ce point de l'objet aura donc pour expression [Bittoun,
  1984], [Brent, 1986] :
  (1.28) Ce champ effectif forme avec l'axe ox un angle j3 donné à chaque point de l'objet par:
  32 f3(X,Y,Z)
  = -arctan(YB~
  -OJ
  RF
  yBt
  J
  (1.29)
  A la fin de l'impulsion, tout se passe donc comme si l'aimantation avait tourné autour de la direction du champ effectif d'un certain angle 8 par rapport à sa direction initiale [Liang]. Cet angle est donné par :
  (1.30)
  avec I10J = yB~ -
  OJ RF •
  La relation 1.25 nous permet de simuler n'importe quelle rotation du vecteur de l'aimantation locale en appliquant l'opérateur
  RRF'
  Pour décrire l'évolution de cette aimantation pendant la période de la précession libre, nous devons ici faire appel à l'opérateur de relaxation
  Rrelax'
  A chaque instant t ' pendant cette
  période, le mouvement de l'aimantation de chaque voxel est le résultat d'un mouvement de rotation et d'un phénomène de relaxation. La rotation s'effectue autour de l'axe oz et elle est définie par l'angle:
  ('
  f
  O(X,Y,Z,t ' ) = I1OJ(X,y,Z,R)dR
  (1.31 )
  o
  Quant au phénomène de relaxation, il fait apparaître les deux fonctions d'exponentielles amorties:
  -l'
  El (t ') = e 11
  (1.32) (1.33)
  L'opérateur de relaxation s'exprime ainsi par la matrice suivante [Bittoun, 1984], [Brent, 1986] :
  33 -E 2 sin(B)
  E 2 cos(B)
  (1.34)
  o Le vecteur aimantation à l'instant t' sera de la forme:
  m(X,Y,Z,t') = Rre1ax(X,Y,Z,t')* m(X,Y,Z,t + r) +m o(X,Y,Z)(1-E 1 )*k
  (1.35)
  Le signal RMN à l'instant t'est obtenu en déterminant le vecteur d'aimantation de chaque voxel de l'objet grâce à cette dernière équation et en effectuant la sommation de toutes les composantes transversales
  A noter que les opérateurs d'impulsion RF et de relaxation ont été présentés dans de nombreuses publications ([Bittoun, 1984], [Brent, 1986], [Brenner, 1997]) et ne seront pas détaillés davantage dans ce document. Il s'agit uniquement pour nous de les utiliser dans la conception du simulateur de signaux RMN.
  V.3.
  Résultats de simulation: exemples de signaux RMN
  V.3.1. Simulation d'une expérience avec une impulsion rectangulaire: L'expérience de base en RMN consiste à appliquer une impulsion rectangulaire RF de durée r et d'angle de basculement a, et à observer la précession libre de l'aimantation globale (signal FID) pendant un temps T. Le type du signal observé à savoir une exponentielle amortie ne permet pas de révéler la validité du simulateur mis en place. Pour pouvoir valider notre simulateur, nous nous sommes basés sur une expérience particulière décrite dans la littérature. Dans cette expérience, un gradient est appliqué pendant l'excitation et le signal est observé en conservant ce même gradient pendant la période de la précession libre. Le signal obtenu approche la forme de l'excitation et sa durée effective (temps de relaxation effectif) étant approximativement égale à la durée de l'impulsion [Mansfield 1989].
  34 Nous avons donc considéré, dans notre simulation, un objet ID homogène de dimension 25cm. Cet objet est discrétisé en petits éléments auxquels on attribue des valeurs identiques d'aimantation mo (densités identiques de protons). Les temps de relaxation longitudinale et transversale sont respectivement
  T = 100ms
  et T2
  = 10ms
  pour tous les
  points de l'objet. La valeur du champ magnétique Bo est fixée à 0.1 T et celle du gradient selon la direction de l'objet à 2.5mT/m. Nous supposons que les inhomogénéités de ce champ sont négligeables. Le signal RMN que nous obtenons en réponse à une impulsion rectangulaire 90° de durée Ims est donné dans la figure 1. 7. La forme et la durée effective de ce signal comparées à celles de l'excitation appliquée confirment le potentiel de notre simulateur à reproduire les observations décrites ci-dessus [Asfour-JMR, 2000] 
  1.2 . -- - - - - - - -, -- - - - - - - - , - - - - - - - -- ,
  Excitation Signal FID
  0.6
  0.4
  0.2
  0 1 - - - - - - - - ' .*  
  -0.2 ' - - - --
  o
  - --
  --"-----------'-------------'
  500
  1000
  1500
  Numéro de l'échantillon n Figure 1.7 Signal simulé de résonance (traits continus) dans le cas d'une impulsion rectangulaire (en pointillé) 90° de Ims de durée et en présence d'un gradient de 2.5mT/m pendant l'excitation et pendant l'acquisition du signal.
  35 V.3.2. Simulation d'une expérience avec une impulsion de forme quelconque: Tel qu'il a été décrit dans le paragraphe V.2, l'opérateur RRF de l'impulsion ne prend pas en compte les variations dans le temps de l'amplitude
  Bt (t)
  du champ radiofréquence
  BI'
  De par sa définition, cet opérateur s'applique uniquement au cas d'une impulsion rectangulaire. La simulation des impulsions de forme quelconque consiste à décomposer l'impulsion en une suite d'impulsions rectangulaires d'amplitudes différentes [Asfour-JMR, 2000], [Bittoun, 1984], [Brenner, 1997], [Olsson, 1995]. L'action de chaque impulsion rectangulaire sera alors décrite par un opérateur matriciel. Nous disposons donc d'un ensemble d'opérateurs R~F à appliquer successivement à l'aimantation locale. Cette dernière après l'impulsion i
  s'écrit:
  -i m
  où
  m
  _ i 1
  = RiRF
  -
  (I.3 6)
  • mi-!
  est le vecteur d'aimantation après l'application de l'impulsion rectangulaire i-l.
  L' application de cet ensemble d'opérateurs est donc équivalente à l'application d'une seule impulsion modulée en amplitude. Typiquement, nous avons décomposé l'impulsion modulée en 1024 impulsions rectangulaires. L'intérêt de cette technique de décomposition réside dans sa capacité d'obtenir une solution numérique de l'équation de Bloch dans le cas général. Plusieurs formes de modulation d'amplitude ont été testées. Celles-ci comprennent notamment les impulsions de forme gaussienne, sinus cardinal (sine) ou encore les impulsions modulées en bruit aléatoire. La figure 1.8 montre le signal RMN correspondant à l'objet ID en réponse à une impulsion sine 90° de durée de Ims. Un gradient de 2.5mT/m a été appliqué durant l'excitation et durant la période de la précession libre. Sur cette figure, nous constatons, encore une fois, que la forme du signal RMN approche la forme de l'excitation [Asfour-JMR, 2000]. Le temps de relaxation effectif T2• étant du même ordre de grandeur que la durée du lobe principal du sine.
  36 Excitation RF 0.8
  SignalRMN
  0.6
  0.4
  0.2
  o :". :" , , " , , ' ,
  .".
  ;:
  ::
  -0.2
  o
  100
  200
  300
  400
  500
  600
  700
  800
  900
  1000
  Numéro de l'échantillon n Figure 1.8 Signal simulé (trait continus) de résonance en réponse à une impulsion sinc 90° (en pointillé) d'une durée de Ims en présence d'un gradient de 2.5mT/m pendant l'excitation et pendant l'acquisition du signal.
  Le travail présenté au cours de ce chapitre a conduit au développement et à la validation d'un simulateur de signaux RMN. De par les principes algorithmiques retenus, ce simulateur permet d'implanter, non seulement une expérience RMN de base, mais toute séquence d'imagerie. Il permet donc la réalisation d'images simulées pour étudier l'influence des paramètres de l'expérience sur la qualité de l'image. Ce simulateur se veut un réel outil offrant un complément d'étude à l'expérience pour toutes équipes de recherche traitant de la RMN. Pour la souplesse d'utilisation de cet outil, nous envisageons le développement d'une interface graphique appropriée.
  37 Chapitre Il
  Description et caractéristiques de l'imageur développé pour l'IRM dédiée à bas champ
  38 1.
  Introduction Dans ce chapitre, nous décrivons la conception et le développement du nouvel imageur
  numérique d'IRM que nous avons élaboré durant cette thèse. Cette description est complétée par l'évaluation des performances des différents constituants de la chaîne d'acquisition des signaux. Nous rappelons qu'une installation typique d'imagerie par résonance magnétique comporte de loin deux éléments majeurs: l'aimant créant le champ magnétique statique
  BQ'
  et l'électronique du traitement et l'instrumentation associée. La création du champ
  magnétique ne fait pas l'objet de notre travail. Il s'agit uniquement pour nous de réaliser une chaîne électronique complète permettant l'émission et la réception de signaux RF, le traitement numérique des signaux RMN et la reconstruction de l'image. Cette électronique est architecturée autour d'une carte d'acquisition et de traitement du signal temps réel (DSP). L' originalité de l' instrumentation développée réside dans la spécification et l'intégration d'un récepteur et d'un synthétiseur de fréquence entièrement numériques. Ce choix instrumental trouve sa justification dans notre volonté de développer un imageur travaillant à bas champ et dédié à l'imagerie de petits organes. Outre son coût peu élevé [Gries, 1987] qui permet son utilisation dans le contrôle agroalimentaire [Constantinesco, 1997], le bas champ est particulièrement avantageux pour la technique dite d'IRM interventionnelle [Ravier, 1996]. Dans ce contexte du bas champ où le rapport signal sur 7
  bruit (RSB) est très faible (RSB oc BQ4), le choix d'une solution numérique permet, contrairement aux systèmes analogiques commercialisés, une amélioration du RSB. Cette solution présente également l'avantage d'être moins coûteuse et moins encombrante [Raoof, 2002].
  Il.
  Description de l'imageur numérique La figure 11.1 présente l'architecture autour d'un DSP du système d'acquisition et
  de traitement numérique des signaux RMN. A l'exception de l'électro-aimant créant le champ magnétique statique
  BQ
  et les bobines de gradients associées, nous avons développé
  l'imageur dans son intégralité. Ceci comprend donc la conception de l'antenne et du duplexeur pour l'émission des impulsions d'excitation et pour la détection de signaux
  39 RMN; de l'amplificateur de puissance RF; du synthétiseur numérique et du récepteur numérique.
  Synthétiseur numérique
  puissance
  4.783 MHz
  40MHz BUS
  Sorties CNAs CarteDSP }
  PC
  BUS
  Duplexeur
  Electroaimant et bobines de gradients
  -
  ver s les ampliOcateur s des gradients Signal RMN
  Port série duDSP
  Récepteur numérique 40MHz
  Figure II.1
  Architecture générale du système numérique pour l'IRM
  Préalablement au développement de cet imageur, nous avons, bien entendu, conduit une étape de recherche instrumentale qui a relevé que les systèmes analogiques disponibles sur le marché ne peuvent pas répondre à nos besoins spécifiques en terme du RSB de récepteur, de la résolution en fréquence et en phase du synthétiseur, du coût, et de l'encombrement.
  Le principe général de l'imageur numérique illustré par la figure II.I s'articule autour la génération, par le synthétiseur numérique, de l'impulsion d'excitation à la fréquence de travail fixée dans notre cas à 4.783MHz (ce qui correspond à une valeur de champ d'environ 0.11 T). Cette impulsion, pouvant être modulée en amplitude, en fréquence ou en phase, est ensuite amplifiée grâce à un amplificateur RF linéaire de puissance. Une fois amplifiée, elle est transmise à l'antenne via le duplexeur qui isole le récepteur pendant cette période d'excitation. L'antenne envoie, après la fin de l'impulsion,
  40 le signal RMN au récepteur numérique via ce même duplexeur. Le récepteur numérique assure la démodulation numérique du signal RMN reçu dont les échantillons dans la bande de base sont transmis au DSP via son port série synchrone à une vitesse de transmission de 16Mbits/seconde. Pour assurer la synchronisation de tous les évènements temporels conduisant à l'obtention d'un signal ou d'une image, le système est cadencé à la fréquence de l'horloge externe du DSP (40MHz). Cette d'horloge fournit notamment la fréquence de référence du synthétiseur numérique qui, à son tour, délivre un signal d'horloge pour le pilotage du récepteur numérique. La fréquence de cette dernière horloge détermine en réalité la fréquence d'échantillonnage du signal dans le récepteur numérique. Le contrôle des trois gradients lors d'une séquence d'imagerie est assuré par trois convertisseurs numérique-analogique (CNA) du DSP.
  Nous allons dans la suite de ce chapitre présenter les différents éléments de ce dispositif expérimental.
  Il.1.
  L'électroaimant
  C'est un aimant résistif 1 produisant un champ magnétique statique vertical Bo d'environ 0.11T dans un entrefer de 15cm. Il dispose également des bobines, créant des gradients de champ magnétique dans les trois directions de l'espace, XYZ, et des bobines de correction de l'homogénéité du champ principal
  Bo.
  Nous avons effectué des mesures pour évaluer la meilleure homogénéité pouvant être obtenue selon les trois axes X,
  y, Z Les résultats ont montré qu'une homogénéité
  d'environ 20ppm (partie par million) est assurée pour chaque axe sur une distance de
  ± 5em par rapport au centre de l'aimant.
  Il.2.
  L'antenne et le duplexeur
  Dans un système d'IRM, le champ radiofréquence tournant BI (ou impulsion d'excitation) est créé par des enroulements de fils conducteurs, appelés antennes, parcourus par des courants électriques sinusoïdaux. Ces enroulements sont conçus pour 1
  Drusch, France. Mis à notre disposition par le Laboratoire de Biomécanique : CHU de Strasbourg.
  41 créer un vecteur champ BI orthogonal au vecteur champ Bo et identique dans tous les points du volume à imager. Cette antenne est chargée généralement de recevoir le courant électrique traduisant le signal RMN. Il s'agit d'antenne émettrice-réceptrice2 . Lors de la conception d'une antenne pour l'IRM, certains critères fondamentaux doivent être pris en compte pour garantir une émission et une réception optimales. Pour créer un champ BI dont la composante active tourne à la vitesse angulaire CO o
  = 27ifo,
  l'antenne doit être accordée à la fréquence de Larmor imposée par Bo' Cet
  accord est assuré généralement par une capacité C qui forme avec l'antenne un circuit parallèle RLC résonant dont la fréquence de résonance est alors:
  1 co = - o .JLC
  (11. 1)
  où L est l'inductance de l'antenne. Le facteur de qualité Q de ce circuit, qui est donnée par la relation 11.2, traduit la sélectivité de l'antenne en terme de fréquence.
  Q=~
  (11.2)
  R
  où R est la résistance série de l'antenne. Plus ce facteur de qualité est élevé, plus la résonance est aiguë et plus l'intensité du courant parcourant l'antenne est élevé, ce qui est appréciable en particulier pour la réception du signal RMN. Le champ BI crée par l'antenne doit être homogène dans tout le volume à imager afin d'assurer une excitation identique de l'intégralité de l'échantillon à étudier. Il faut ajouter aussi qu'en imagerie, on peut montrer que l'effet des inhomogénéités de BI se traduit par des pertes du signal (perte en sensibilité de l'antenne) [Pimmel, 1990].
  La structure et la géométrie d'une antenne RMN varient énormément d'une application à l'autre. Elles doivent cependant être adaptées à l'échantillon à étudier [Choquet, 1996]. Dans notre application, nous avons choisi la structure d'antenne en 2
  On trouve aussi des systèmes où l'émission et la réception sont assurées par deux antennes distinctes. On se
  limite ici au cas des antennes assurant un double rôle d'émission et de réception.
  42 solénoïde qui est particulièrement adaptée à l'imagerie de volume. Cette structure est parmi les plus simples à réaliser et a été longuement étudiée et développée dans des travaux antérieurs [Choquet, 1996], [Pimmel, 1990]. Un solénoïde est une inductance formée par un conducteur de rayon r enroulé en hélice (n spires) autour d'un cylindre de diamètre d et de longueur 1. L'échantillon est placé à l'intérieur du solénoïde. Les principaux paramètres géométriques (d, l, r et nombre de spires n) de ce type d' antenne sont calculés pour obtenir un facteur de qualité le plus élevé possible tout en ayant un champ
  BI
  homogène et suffisamment intense. On montre qu'il n'existe qu'un
  compromis entre le facteur de qualité et l' homogénéité de
  BI ' Les paramètres de l' antenne
  permettant de réaliser l' optimal sont spécifiés sur le schéma d'antenne présenté dans la figure 11.2.
  1 ~ 0.8d
  d
  1 1
  Figure D.2
  1 :
  1
  1
  -.::+-
  -.:
  :.-
  2r
  S =
  3r à 4r
  1
  1
  Schéma d'une antenne solénoïde de dimensions optimales
  Nous avons réalisé une antenne solénoïde de dimensions avoisinant celles données sur cette figure. Les caractéristiques principales de cette antenne sont résumées dans le tableau IV.I. La fréquence de résonance et le facteur Q de qualité ont été mesurés par la méthode de la bande passante.
  43 Diamètre
  Longueur
  rayon du
  Nombre de
  Fréquence de
  Facteur de
  d
  1
  conducteur
  spIres n
  résonance
  qualité Q
  mesurée
  mesuré
  4,783MHz
  190
  r
  4cm
  3.4cm
  1.5mm
  Tableau n.1
  7
  Caractéristiques de l'antenne solénoïde.
  Les signaux RF d'émission et de réception sont couplés à l'antenne par une boucle, dite d'accord, par l'intermédiaire du duplexeur comme le montre la figure 11.3.
  r*- * *-I~* _*~: ;J:':~::::~_. . .*_-_*_- _*_*_. .*i
  !
  '1
  :
  ~
  ' : S2
  duDSP
  ,i
  ! i
  o "~ Antenne et capacité d'accord
  !,
  Commande TIL du DSP
  t
  Câble/oaxi al
  "4
  SI c******* ****i / - - - - - - j!r----"
  0;
  Ll."-:"mJ
  i 1
  i ~:\!:.';' 1
  , Boucle !. d'accord 1:
  : i:
  L_ __• ____ __  _ ____    _ _.   _ _ __  ____  __  __ ___ 1
  Le duplexeur
  Figure II.3
  Schéma électrique équivalent de l'antenne et du duplexeur.
  La boucle d'accord schématisée sur cette figure par son inductance permet de modifier la fréquence de résonance et l'impédance de l'antenne et ce en raison des modifications de l'inductance mutuelle pouvant être générées par les déplacements de la boucle par rapport à l'antenne. Ceci constitue un moyen simple permettant à l'antenne d'être à nouveau accordée et réadaptée en impédance. Ces deux caractéristiques étant modifiées par l'introduction dans l'antenne d'un échantillon. 44 L'antenne étant utilisée en émission et en réception, le duplexeur a pour rôle d'assurer la séparation du circuit d'émission de celui de la réception. Pendant la phase d'émission, celui-ci est chargé de transmettre l'impulsion d'excitation à l'antenne tout en protégeant l'entrée du récepteur numérique de ces mêmes impulsions. En revanche, durant la période de réception, le duplexeur permet la transmission du signal RMN en provenance de l'antenne vers le récepteur tout en assurant l'isolation efficace du circuit de réception de celui de l'émission. Cette séparation évite une détection de bruit provenant de ce dernier. Le duplexeur est, en fait, composé principalement de deux interrupteurs RF, de deux diodes de signal montées en tête-bêche et d' un câble coaxial de longueur Â" soit 4 environ 16m pour la fréquence 4.783MHz. Les deux interrupteurs RF sont commandés par des signaux logiques compatibles TTL générés par le DSP. En émission, les signaux de commande sont tels que les interrupteurs 1 et 2 sont respectivement fermé et ouvert. Les diodes sont conductrices le l'excitation RF est transmis à la boucle d'accord. Au niveau du récepteur, l' extrémité du câble coaxial
  Â,
  4
  se trouve en court-circuit assuré par l'interrupteur 2 alors ouvert.
  L'impédance ramenée au niveau de la boucle d'accord est infinie et la réception se trouve par conséquent verrouillée. Aucun signal ne parvient donc au récepteur. A la fin de l'impulsion, i.e. en réception, le rôle de deux interrupteur s'inverse et les diodes sont donc bloquées. Le signal RMN détecté par la boucle d'accord parcourt donc le seul chemin possible vers le récepteur via le câble quart d'onde.
  Il.3.
  L'amplificateur RF de puissance
  Il est chargé d'effectuer l'amplification nécessaire pour la production d'un champ BI suffisamment intense. Les caractéristiques les plus relevantes de l'amplificateur RF de puissance pour l'IRM sont la linéarité, la sélectivité autour de la fréquence d'intérêt et la puissance maximale disponible. Cette dernière est d'environ 200W pour notre imageur. Les détails de conception et de réalisation de cet amplificateur ne seront pas exposés ici. Nous présentons brièvement ses principales composantes et nous donnons également ses caractéristiques de performance liées à notre objectif.
  45 L'amplificateur que nous avons développé est fonné de l'association de deux étages d'amplification réalisés sur deux circuits imprimés indépendants. Le premier est un préamplificateur large bande de faible puissance. L'amplification de puissance proprement dite est réalisée par un deuxième étage sélectif autour de la fréquence du travail. Il.3.1. Le préamplificateur
  Il est constitué par la mise en cascade de trois étages d'amplification à base de transistors de puissance bipolaires (figure II.4). Le premier étant un montage émetteur commun fonctionnant en classe A et réalisant aussi un filtre passe-haut dont la fréquence de coupure voisine le 3MHz. Le deuxième étage, fonctionnant également en classe A, a la même structure que le premier. Il possède à son entrée un filtre passe-bas de fréquence de coupure d'environ 20MHz. Ces deux étages fixent ainsi la bande passante du préamplificateur à 3-20MHz et pennettent d'amplifier le signal d'entrée (généralement un petit signal) à un niveau nécessaire à l'attaque du troisième étage. Celui-ci est composé d'un amplificateur travaillant en classe B fonné par deux transistors montés en push-pull pour assurer l'amplification finale du signal.
  Vers l'étage sélectif
  ~-----------
  Y ----------_/
  Deux étages classe A
  ~-------V -------~ Etage Push Pull classe B
  Figure II.4
  Schéma bloc du préamplificateur.
  46 Les caractéristique fréquentielles du préamplificateur sont présentées dans la figure
  n.5 en terme de son gain en tension Av. Dans la bande de fréquences 3-20MHz, un gain en tension d'environ 24dB est obtenu.
  26~---.----.----.----.----.----.----.----~----.----.---.
  24
  22
  co ~Q) 20
  ~
  t5 18
  OUi
  c
  2 c
  Q)
  o~ 16 (!)
  14
  12
  10 ~--~----~--~----~--~----~--~----~----~--~--~
  o
  2
  468
  10
  12
  14
  16
  18
  20
  22
  Fréquence en MHz
  Figure II.5
  Gain en tension du préamplificateur en fonction de la fréquence.
  Il.3.2. L'étage sélectif Il est constitué d'un amplificateur de puissance fonctionnant en classe C autour de la fréquence d'intérêt. Nous avons choisi de réaliser ce type d'amplificateur en raison de l'adéquation de ses caractéristiques à nos besoins. D'une manière générale, un amplificateur de classe C a un rendement supérieur à celui de la classe B. Toutefois, pour amplifier un signal sinusoïdal, il faut l'accorder à la fréquence de celui-ci. L'amplificateur classe C est dit sélectif et réalise l'amplification dans une bande étroite autour d'une fréquence centrale.
  47 Le schéma bloc de cet étage d'amplification est représenté sur la figure 11.6. Le circuit résonant placé à la sortie du transistor permet de sélectionner les fréquences dans une bande étroite autour de sa fréquence de résonance 5.1MHz. Ce circuit offre la possibilité de pouvoir accorder l'amplificateur à n'importe quelle autre fréquence de résonance dans la bande passante du préamplificateur. Ceci nous permet d'adapter l'amplificateur de puissance pour d'autres fréquences de travail (donc pour d'autres valeurs du champ magnétique statique) dans la bande 3-20MHz. Actuellement, une autre version de cet amplificateur est en cours de développement pour Un imageur travaillant à 15MHz.
  Sortie d u préamplific ateur
  
  ~
  
  Circuit résonant
  
  
  Vers l'antenne
  '--------------~-------------~ Etage sélectif classe C Figurell.6
  Schéma bloc de l'étage sélectif
  1/.3.3. Caractéristiques mesurées de l'amplificateur total Les caractéristiques fréquentielles de l'amplificateur complet sont représentées sur la figure II. 7. La fréquence de résonance est voisine de 5.1 MHz. Pour cette fréquence, le gain total en tension est d'environ 43dB, soit 19dB de gain pour l'étage sélectif.
  48 45 40 35 30 []J
  "0 ;:: QI
  ~;:: .~
  25 20
  ;::
  .2l ;:: QI
  ;::
  15
  cii
  (9
  10 5 0 -5
  0
  Figurell.7
  2
  4
  6
  8 10 12 Fréquence en MHz
  14
  16
  18
  20
  Gain en tension de l'amplificateur RF en fonction de la fréquence
  La figure IL8 présente les variations de la tension de sortie en fonction de celle de l'entrée qui montrent la linéarité de l'amplificateur dans la bande passante. Sur cette même figure est donnée aussi la puissance crête à crête en sortie en fonction de celle de entrée. Une puissance maximale crête à crête d'environ 180W peut être donc obtenue.
  49 100.------,,------.-------,------,-------,-------,------, 90
  ~
  80
  Q)
  ~
  70
  fi)
  c: Q)
  8
  'iii
  '5
  60
  ~ Q)
  'C
  2 'Q) t -lU Q)
  ~ (,) Q)
  'C
  :ê
  a. ~
  20
  10
  0.1
  0.2 0.3 0.4 0.5 Amplitude crête à crête de l'impulsion en entreé (Volt)
  0.6
  0.7
  160.-----.------.-----.------.-----,,-----.------.-----.-----~
  160
  ~ Q)
  'E
  140
  fi)
  c: Q)
  8
  120
  'iii
  '5 a.
  ~ 100 Q)
  u
  Q)
  ~
  60
  (,)
  'aI
  2
  '~ 60 Q) (,)
  c: al
  :ll
  '5
  40
  c 20
  2 3 4 5 6 7 Puissance crête à crête de l'impulsion en entreé (mW)
  8
  9
  Figure II.8 Caractéristiques entrée-sortie de l'amplificateur RF. En haut: l'amplitude crête à crête de l'impulsion en sortie en fonction de l'amplitude crête à crête d'impulsion en entrée; En bas: la puissance crête à crête de l'impulsion en sortie en fonction de celle en entrée (puissances mesurées sur une charge de 50n.
  50 Il.4.
  La Carte DSP (Digital Signal Processor)
  Elle constitue le coeur du système d'imagerie. Cette carte 3 d'acquisition, de contrôle et de traitement est architecturée autour du processeur de traitement du signal DSP (Digital Signal Processor) ADSP21061 SHARC (Super Harvard Arcitecture Computer) d'Analog Devices. C'est une carte de type «Plug and Play» interfacée avec le PC via le bus PCI de celui-ci. Comparé à un processeur conventionnel, l'utilisation de DSPs dans le système d'IRM offre plus de possibilités et de flexibilité en traitement du signal et en reconstruction d'images [Morgan, 1999]. Les processeurs ADSP21061 sont des processeurs 32 bits travaillant en virgule flottante. Ils sont optimisés pour le traitement numérique du signal en temps réel et largement utilisés dans les applications nécessitant une grande vitesse de calcul. Leur réputation vient principalement de leur capacité d'effectuer les calculs en virgule flottante. Ils répondent ainsi à des critères de qualités plus élevés que leurs équivalents en virgule fixe. L'horloge interne (160MHz) du processeur est générée à partir d' une horloge externe sur la carte. Cette dernière horloge de fréquence égale à 40MHz constitue également l ' horloge de référence du synthétiseur numérique (voir fi gures il. 1). Le processeur possède une mémoire interne de 4Mbits (4mégabits) et dispose de plusieurs protocoles rapides et efficaces pour la communication et le transfert de données avec une large variété de périphériques (mémoire externe, PC, autres DSP, périphériques de traitement numérique du signal.  ). En particulier, nous utilisons deux ports série synchrones pour le transfert de données numériques (échantillons du signal RMN) du récepteur vers la mémoire interne du DSP. Des bus parallèles nous permettent également de charger, à partir du DSP, les paramètres du récepteur numérique et du synthétiseur numérique (voir figures il.1, Il.9 et il.11). En plus du processeur, la carte, intègre un convertisseur analogique-numérique (CAN) et trois convertisseurs numérique-analogique (CNA). Le CAN possède 4 entrées qui peuvent être converties les unes après les autres sur 12 bits à une fréquence d'échantillonnage paramétrable entre 153Hz et 29kHz Les trois CNAs transforment des données à virgule flottante en des données à virgule fixe avant de procéder à la conversion en signal analogique dont le niveau est réglable dans la gamme ± 10 Volts. Le bruit résiduel à la sortie est très inférieur à 1
  51 mVolt. Ces trois CNAs sont chargés de contrôler les trois gradients lors d'une séquence d'imagerie (figure II.1).
  La programmation de la carte peut se faire en langage assembleur ou en langage C. Pour optimiser le temps du calcul et l'utilisation des mémoires, nous avons choisi de développer tous nos programmes en assembleur.
  Pour plus de détails concernant l'architecture du processeUr ADSP21061 et les possibilités offertes, le lecteur pourra consulter le document de synthèse donné en annexe.
  Il.5.
  Le synthétiseur numérique
  Il.5.1. Intérêt du synthétiseur numérique Dans un instrument d'IRM, les tâches que doit accomplir le synthétiseur de fréquence se multiplient en raison du développement rapide de nouvelles techniques (séquences) d'imagerie. Ces techniques, qui sont de plus en plus sophistiquées, placent des exigences supplémentaires très strictes sur le synthétiseur de fréquence. Celles-ci concernent essentiellement la phase et la fréquence des impulsions générées par le synthétiseur. En effet, quelques méthodes d'imagerie rapides, comme certaines variantes de la séquence d'équilibre dynamique, nécessitent l'incrémentation de la phase des impulsions successives par des quantités programmables [Liang, 1999]. D'autres séquences font appel à une modulation, par un signal modulant prédéfini, de la phase de l'impulsion d'excitation. Plus classique, nous pouvons citer l'imagerie multi-coupes qui impose une modification de la fréquence des impulsions successives et l'imagerie dite décalage du champ d'imagerie 4 (offset field ofview imaging) dans laquelle la fréquence du signal utilisé pour démoduler le signal RMN reçu, doit être décalée par rapport à la fréquence de Larmor. Nous pouvons aussi rajouter certaines techniques d'imagerie qui utilisent une modulation de fréquence de l'impulsion d'excitation [Stormont, 1990]. Plus généralement, le point commun entre toutes les méthodes d'imagerie consiste dans la nécessité de maintenir une relation connue entre la phase du signal d'excitation et
  3 4
  SRT (Signal Recherche & Technologie) : France Voir Chapitre III pour la définition du champ d'imagerie.
  52 celle du signal utilisé pour la démodulation des signaux RMN et ceci pour faciliter, après acquisition, l'application, si cela est nécessaire, d'une correction de phases. Les synthétiseurs actuels peuvent difficilement satisfaire l'ensemble de ces exigences et s'ils le font, ce n'est qu'avec complexité et peu de souplesse. Ces contraintes nous ont conduit à choisir une solution technologique consistant dans l'utilisation d'un synthétiseur numérique programmable avec une grande flexibilité et aisance. Le synthétiseur numérique AD9852 que nous avons donc intégré dans notre système d'imagerie constitue une grande nouveauté dans le développement instrumental de l'IRM. Il peut générer le signal de l'excitation avec une fréquence et une phase précises et quelconques. Ce signal peut être facilement modulable en amplitude, en fréquence et en phase.
  Il.5.2. Le fonctionnement du synthétiseur AD9852 Le synthétiseur numérique AD9852 d'Analog Deviees est un circuit hautement intégré qui utilise la technologie avancée de la synthèse directe de fréquence encore appelée DDS (Direct Digital Synthesizer). Cette technologie innovatrice et ultra rapide ne fait pas appel à l'usage de boules à verrouillage de phase (PLL) systématiquement employées dans les synthétiseurs de fréquence. En effet, les synthétiseurs de fréquence basés sur l'utilisation d'un PLL présentent les inconvénients d'instabilité face au changement rapide de fréquence. Ils souffrent également de la présence d'un bruit de phase non négligeable [Raoof, 2002]. La technologie DDS est basée sur le fait qu'on peut définir une sinusoïde, dont on connaît la fréquence et la phase, en spécifiant une série de valeurs (échantillons) prises à des intervalles égaux définis par la fréquence d'échantillonnage. Cette série de valeurs est envoyée à un convertisseur numérique-analogique qui fournira donc le signal analogique.
  La figure 11.9 montre le diagramme simplifié du synthétiseur AD9852 utilisé.
  5
  M. Decorps : notes de cours DEA Méthodes Physiques Expérimentale, 1998.
  53 DDS
  /
  /
  -'<-.,>' /
  I2-bits
  CNA 1 I2-bits
  -
  /
  Sortie du synthé,!iseur 
  a I2-bits
  / /
  //
  .4l
  Jl
  32-bits
  r+
  / / I4-bits ./
  Fréquence
  amplitude
  Phase
  CNA2 I2-bits
  
  V lEntrées .-_ _ _--tt-jnalogiques
  L.~~.
  4 Registres de programmation
  Sortie analogique
  _ . /
  V 1
  
  Sortie horloge
  comparateur
  /
  Horloge de référence (40MHz) ~
  V Liaison parallèle 8 bits Chargement des paramètres (fréquence, phase, amplitude par le DSP
  ** 
  ____~.
  Figure II.9
  DSP **••••••. :
  Diagramme simplifié du synthétiseur numérique AD9852
  Le DDS, associé à un convertisseur numérique-analogique (CNA 1) interne de 12 bits ultra rapide (300MHz), fonne un circuit synthétiseur numériquement programmable. Référencé à une source d'horloge stable, ce synthétiseur produit une sortie très stable en tenne de fréquence, de phase et d'amplitude. Fréquemment utilisé dans les applications de télécommunications, ses caractéristiques remarquables sont des atouts majeurs pour son utilisation dans un système d'IRM car elles répondent aux critères et aux exigences mentionnées au paragraphe précédent. L'horloge de référence du synthétiseur est fournie par le DSP (40MHz). Le synthétiseur génère ensuite sa propre horloge interne de fonctionnement. La fréquence de cette dernière est programmable et peut aller jusqu'à 300MHz. Cette fréquence, qui constitue la fréquence d'échantillonnage du signal numérique à la sortie du DDS, est
  54 actuellement fixée 200MHz. Ainsi, la fréquence maximale du signal généré par le synthétiseur est alors de 100MHz. La fréquence, la phase et l'amplitude du signal synthétisé sont directement programmables par le DSP avec la possibilité de mise à jour rapide de la fréquence (100 millions de nouvelles fréquences par seconde). La résolution en fréquence étant de 48 bits, ce qui correspond à 1microHertz pour une fréquence d'horloge interne de 300MHz. Cependant, seuls 32 bits sur 48 bits sont actuellement utilisés pour conserver la même résolution fréquentielle entre le signal d'excitation et le signal utilisé pour la démodulation dans le récepteur et qui possède une résolution de 32 bits (voir paragraphe II.6.2.2). Quant à la phase du signal synthétisé, elle est programmable avec une résolution de 14 bits, soit donc 0.02°. A titre de comparaison, nous signalons que les résolutions en fréquence et en phase sont respectivement de O.IHz et de 0.225° pour un synthétiseur de fréquence utilisé actuellement sur une console conventionnelle d'IRM 6
  Le signal numérique synthétisé par le DDS est envoyée sur le CNA 1 qui délivre le signal analogique (une sinusoïde).
  La disponibilité d'un deuxième convertisseur CNA 2 et d'un comparateur ultra rapide (figure 11.9) est une particularité de grande importance pour notre système d'imagerie dans le contrôle de la phase absolue, i.e. le contrôle de la référence de phase des impulsions lors d'une séquence d'imagerie. En effet, la phase de la première impulsion du cycle de la séquence définit la phase de référence. Cette phase doit être conservée d'un cycle à l'autre durant toute la séquence. Pour ce faire, cette phase de référence peut être définie par un nombre programmable par le DSP et envoyé au CNA 2 qui donne un niveau de tension fixe à sa sortie. L'envoi d'une impulsion ayant la même phase de référence se fait en comparant, grâce au comparateur, la phase actuelle de la sortie du synthétiseur avec la tension de référence. Le résultat de comparaison est interprété par le DSP qui sera chargé d' ajuster la phase de la sortie du synthétiseur pour avoir la valeur de la phase de référence.
  6
  Données techniques du fabriquant S.M.I.S.
  55 Il.6.
  Le récepteur numérique
  Comme nous l'avons mentionné au premier chapitre (paragraphe II2.3), le signal RMN est généralement détecté en quadrature. La performance de cette détection est essentiellement caractérisée par le rapport signal sur bruit (RSB) du récepteur. Celui-ci constitue donc l'élément clé qui doit permettre d'extraire le signal utile avec un maximum duRSB.
  Il.6.1. Rappel : récepteur analogique La structure typique d'un récepteur utilisé dans un instrument conventionnel d ' IRM est représentée sur la figure II. 10
  CL
  1
  CAN
  Filtre passe b as
  Signal RlvIN
  (t ~ Antenne
  0
  CIl
  Synthétiseur local de fréquence
  '"0
  Préamplificateur Déphaseur
  CL
  Q
  CAN
  Filtre passe bas
  Figure II.10
  Schéma général d'un récepteur analogique.
  Le signal en provenance de l'antenne est amplifié par un préamplificateur de gain fixe (40-50dB). Ce signal est dirigé sur deux voies et multiplié par les deux signaux (sinus
  56 et cosinus) fournis par le synthétiseur local de fréquence. Les sorties des deux voies sont filtrées par un filtre passe-bas pour donner les deux signaux basses fréquences dans la bande de base. Ces deux signaux, correspondant à la partie réelle l et la partie imaginaire Q du signal RMN, sont appliqués respectivement à un convertisseur analogique numérique fonctionnant typiquement à une fréquence d'échantillonnage de quelques 10kHz. Les échantillons en format numérique sont ensuite transmis vers une carte d'acquisition et de traitement DSP pour être exploités dans le calcul du spectre ou la reconstruction de l'image par exemple. Peu de traitement numérique du signal est effectué dans ce type de récepteur Le récepteur analogique a des limitations liées à l'utilisation de composantes analogiques, i.e. l'amplificateur, le synthétiseur local, le déphaseur (7r ), le multiplieur et 2 les filtres. Toute composante analogique est, comme nous le savons, une source du bruit thermique qui dégrade significativement le RSB . Les deux voies de la chaîne analogique, supposées parfaitement identiques, ne peuvent pas l'être en pratique en raison des imperfections des composants analogiques. Par conséquent il peut y avoir un niveau continu sur l'une et/ou l'autre voie, une différence de gain entre elles ou un défaut de déphasage. Ceci est aggravé dans le temps par la sensibilité des composants analogiques aux variations de la température et au vieillissement. L'ensemble de ces défauts, qui ne sont pas facilement maîtrisables, conduit à l'apparition d'artéfacts, appelés artéfacts de quadrature, sur l'image. La détection en quadrature peut être grandement améliorée par l'utilisation de techniques de traitement numérique du signal. Il devient aujourd'hui très attractif d'accomplir le maximum des étapes de la détection en quadrature dans le domaine numérique. On parlera donc de récepteurs numériques (radio logiciel). En d'autres termes, le convertisseur analogique-numérique peut être déplacé en amont de la chaîne de réception au plus près de l'antenne. Cela est rendu aujourd'hui possible grâce aux avancées technologiques dans le développement des convertisseurs analogique-numérique en terme de résolution (nombre de bits) et de fréquence d'échantillonnage [Walden, 1999]. En effet, le développement accéléré de systèmes numériques de traitement du signal haute fréquence est le fruit du marché en pleine expansion de la téléphonie mobile. La demande est tellement forte que le rapport qualité/prix des systèmes numériques ne cesse d'augmenter.
  57 Il.6.2. Intégration d'un récepteur numérique Il.6.2.1.
  Description générale
  Nous avons donc intégré dans notre système d'imagerie un récepteur numérique (figure II.11) typiquement employé dans les applications de télécommunication. Ce récepteur présente plusieurs avantages par rapport au récepteur analogique. Il permet un échantillonnage direct des signaux pour améliorer les performances de la réception et réduire le coût et le bruit de l'électronique. Après échantillonnage, la démodulation et le traitement du signal se font numériquement avec plus de flexibilité et de précision qu'en traitement analogique. Tous les problèmes rencontrés avec le récepteur analogique sont donc résolus ou au moins minimisés grâce à la réception numérique.
  DDC
  1
  3 Nitres de décimation Contrôle automatique du gain
  entrée
  Synthétiseur numérique local
  CAN
  Signal RMN
  Fréquence __~d~'oc~'~h=an=t=ill~o=nn~a~e__~________~~
  Q * 3 Nitres de décimation *** *
  40MHz
  Port série ------tBUS
  DSP Figure II.11
  Schéma du récepteur numérique
  Le récepteur accepte des signaux analogiques et fournit les deux sIgnaux en quadrature (I et Q) en format numérique dans la bande de base.
  58 Le seul élément analogique du récepteur est l'amplificateur à gain numériquement programmable. Ainsi, le signal RMN, encore analogique, est amplifié avant d'être numérisé grâce à un convertisseur analogique-numérique (CAN) 12-bits dont la fréquence d'échantillonnage, numériquement programmable, peut aller jusqu'à 50MHz. Le signal numérique est ensuite démodulé numériquement grâce au démodulateur DDC (pour Digital Down Converter). Dans celui-ci, le signal est multiplié numériquement par les deux signaux de multiplication fournis par un synthétiseur numérique local. Le DDC assure également le filtrage passe-bas grâce à une cascade de trois filtres de décimation de réponse impulsionnelle finie (RIF). Les échantillons numériques de la partie réelle (1) et de la partie imaginaire (Q) du signal sont ensuite transmis au DSP via son port série à une vitesse programmable. Tous les paramètres du récepteur incluant la fréquence d' échantillonnage, la fréquence et la phase du synthétiseur local et les facteurs de décimation des filtres ainsi que leur coefficients sont programmables par l'utilisateur via le DSP.
  Il.6.2.2. Caractéristiques principales du récepteur numérique utilisé o L'amplificateur programmable DVGA (Digital Variable Gain amplifier) C'est un amplificateur à gain numériquement programmable entre -12 et 30dB par pas de 6dB dans une bande passante de 350MHz. Le bruit de l'amplificateur est très faible (2.5nV /.J Hz).
  o Le convertisseur analogique-numérique (CAN) Il est optimisé pour les récepteurs numériques hautes fréquences amSI que pour d'autres applications nécessitant une bonne résolution, une fréquence d'échantillonnage élevée et une large plage dynamique. Un rapport signal sur bruit (RSB) d'environ 67dB est obtenu à la sortie du convertisseur. A noter que l'amplificateur programmable (DVGA) en combinaison avec le DDC peuvent assurer un contrôle automatique du niveau du signal à l'entrée du CAN. Cela permet de
  compresser ou
  d'augmenter la plage dynamique
  du
  signal
  l'échantillonnage. En faisant ainsi, la plage dynamique du CAN est améliorée.
  o Le synthétiseur numérique local
  59
  avant Il est constitué d'un oscillateur numérique local qui fournit les deux signaux de multiplication avec une précision de 32 bits en fréquence, soit une résolution de 0.02Hz pour une fréquence d'horloge de 50MHz, et de 16 bits sur la phase (résolution de 0.005°).
  o Les filtres numériques passe bas de décimation D'une manière générale, un filtre numérique est une alternative au filtre analogique si les performances recherchées dépassent les limites physiques des composants analogiques. Les filtres numériques sont également plus efficaces que leurs *contre parts de filtres analogiques en terme de sélectivité et de linéarité de phase [Chester, 1999]. Les trois filtres numériques du récepteur sont dits de décimation. En plus de son rôle comme filtre passe bas, un filtre de décimation accomplit un sous-échantillonnage du signal. En d'autre terme, si on désigne par
  te la fréquence d'échantillonnage du signal à
  l'entrée d'un tel filtre ayant un rapport de décimation M, le signal à sa sortie se trouve échantillonné à une fréquence égale à
  te. Ceci revient donc à choisir en sortie du filtre un
  M
  échantillon parmi M à son entrée, comme le montre la figure H.12
  ~I
  te
  0
  ~I ~M
  te
  Filtre passe bas
  te M
  •
  1
  te
  le
  2
  2
  Figure II.12
  ~ te
  2M
  illustration du concept de la décimation.
  Ce processus de décimation permet d'augmenter le RSB du convertisseur. En effet, la bande passante d'un système numérique est fixée par la fréquence d'échantillonnage et s'étend de 0 à
  te . Après décimation d'un rapport M, 2
  60
  cette bande se réduit à
  te . Le
  2M rapport signal sur bruit (RSB), qui est inversement proportionnel à la bande passante, est ainsi amélioré. Dans notre récepteur numérique, l'application de cette opération, appelé encore traitement du gain, permet d'améliorer le RSB du convertisseur d'environ 24dB.
  Dans le récepteur numérique, le premier filtre de décimation est du type CIC (pour Cascaded Integrator Comb). TI assure une décimation d'un facteur élevé compris entre 8 et
  2048. Les coefficients de ce filtre sont fixés. Le deuxième filtre possède 21 coefficients programmables par l'utilisateur. Son facteur de décimation est fixé à 2. Ce filtre sert à compenser les ondulations produites par le premier filtre. La réponse impulsionnelle de ce filtre et son gain en dB sont représentés sur la figure ll.13
  35)(10' QG
  . 70
  30
  20 10
  *lO~----!---~ 1 0--"""" ' 5---f2O'------.'
  ~L---~--~10~--~1~'--~20~---.'~
   -
  Figure II.13
  _(>Hz)
  Réponse impulsionnelle et gain du deuxième filtre de décimation.
  Quant au troisième filtre, il a 63 coefficients programmables et un facteur de décimation qui peut prendre une des deux valeurs 2 et 4. Ce filtre est chargé de définir la largueur finale de la bande passante. TI doit donc être adéquatement synthétisé, dans notre cas, de manière à avoir la bande passante désirée pour notre système d'imagerie.
  Il.6.2.3. Evaluation récepteur
  expérimentale
  des
  caractéristiques
  du
  Pour quantifier les performances du récepteur numérique, en terme de RSB, de sensibilité et de réjection de bruit, nous l'avons testé en appliquant à son entrée un signal de niveau connu.
  61 Nous avons, au préalable, fixé la fréquence d'échantillonnage du CAN à 50MHz. La fréquence de démodulation a été fixée à 5MHz. Le facteur de décimation total a été choisi, par les trois filtres, de telle manière que la fréquence d'échantillonnage finale soit de 25kHz dans la bande de base, soit donc un facteur de décimation total de 2000. Nous avons calculé la réponse implusionnelle du troisième filtre de décimation pour avoir une largeur de bande finale de 5.8kHz à - 3dB [Rao of, 2002]. C'est en effet la largeur de bande nécessaire pour les signaux RMN sur notre système. Les coefficients du filtre ainsi calculés sont représentés au même temps que le gain du filtre sur la figure II.14.
  La procédure de test a consisté à appliquer un signal sinusoïdal de fréquence 5MHz et ayant une amplitude réglable entre -50dB et - 130dB. Nous avons constaté que la sensibilité du récepteur était de - 120dB. Nous avons également validé la largueur de bande choisie et en particulier fréquence de coupure du filtre de décimation que nous avons synthétisé. 1 10'
  
  . . . do_
  *'O~---;'~O~-:!20;--~";------;;;40~---:;;50~~"~
  *200~-----C~~~----:':-------C:-----;IO;------;'2;-' _(1Hz)
  Figure II.14 Réponse impulsionnelle et gain du troisième filtre de décimation. La meilleure plage dynamique, en pleine échelle numérique de sortie, a été obtenue pour un signal d'entrée de -96dB, ce qui représente donc un gain du récepteur numérique de 96dB contre environ 80dB pour un récepteur analogique7 • La réjection de bruit hors bande a été évaluée en appliquant au récepteur la somme d'un signal de 5MHz et d'un bruit large bande avec un RSB de OdB. Ce bruit a une largeur de bande de 4MHz de chaque côté en dehors d'une bande de ± 50kHz autour de 5MHz. Les résultats ont montré une réjection totale du bruit hors bande.
  7
  Données techniques du fabriquant S.M.LS
  62 Chapitre III
  Résultats: signaux RMN et premières images
  63 1.
  Introduction Dans ce chapitre, nous allons donner, en premier lieu, quelques exemples des signaux
  les plus caractéristiques en RMN obtenus avec notre système d'imagerie numérique. Dans un deuxième temps, nous présenterons la méthode générale permettant le codage de l'espace et la reconstruction d'une image à partir du signal. La capacité de l'imageur de produire des images sera ensuite illustrée à travers un exemple des premières images obtenues avec quelques essais esquissés basés sur la séquence d'écho de gradient. Nous clôturerons ce chapitre par une brève présentation des travaux en cours de réalisation.
  Il.
  Signaux RMN acquis avec l'imageur numérique Il.1.
  Signal de décroissance d'induction FIC
  La figure III.1 montre un exemple du signal FID démodulé acquis grâce à notre système d'imagerie décrit au chapitre II. Ce signal, obtenu sur un échantillon cylindrique d'eau déminéralisée, est donné par sa partie réelle l et sa partie imaginaire Q qui sont naturellement déphasées de
  7r.
  2
  La fréquence d'échantillonnage est de 10kHz. L'excitation
  étant une impulsion rectangulaire 90° de 2ms.
  64 .
  :
  .
  -1 ',"
  1
  ,
  *1 . . .
  iii"i+i+']i ' !+i'!'+ii'i+iillii+ " +
  oeil* ! 
  Figure m.l
  Il.2.
   l.
  ! 
  Signal de décroissance d'induction FID obtenu avec le système d'imagerie.
  Les signaux d'écho
  Les signaux dits d'écho de spin et d'écho de gradient ont une importance majeure en IRM, car ce sont ces types de signaux qui sont généralement utilisés pour la réalisation de l'image. Nous allons donc les exposer dans ce paragraphe.
  Il.2.1. Signal d'écho de spin Comme nous l'avons vu dans le premier chapitre et l'avons illustré ci-dessus, le signal FID a la forme d'une exponentielle amortie avec une constante de relaxation transversale effective T;. Cette décroissance selon T2• est due, rappelons-le, aux hétérogénéités du champ magnétique statique
  Ba'
  La technique dite d'écho de spin permet d'éliminer à un instant déterminé l'influence des hétérogénéités de ce champ magnétique. Il s'agit, plus précisément, de transformer grâce
  65 à une impulsion RF, un signal d'induction décroissant FID en un signal d'écho [Guinet, 1992]. On se place dans le cas où a été créé, à l'instant t=O, un signal de décroissance d'induction FID, à l'aide d'une impulsion 90° (figure III.2). Si une impulsion 180° est appliquée à un instant T où le signal FID a pratiquement disparu, on constate qu'il apparaît à nouveau un signal d'induction. Ce signal est appelé signal d'écho, i.e. écho du signal FID initial.
  t
  ,
  "\ ,
  Signal FID théorique
  e
  -" "---
  T2 t
  .-_ ' -'",",--Signal FID réel "
  e
  T2•
  t 1 1 1
  0
  .
  : 180°
  ~
  1
  9n
  1'"
  1 1 1 1
  t
  ~I", 1
  T
  ~
  T
  1
  •
  ~
  TE Figure m.2
  Illustration du signal d'écho de spin.
  L'intervalle de temps, ou temps d'écho, noté TE' qUI sépare l'impulsion 90° du maximum du signal d'écho, est égale au double de l'intervalle de temps qui sépare les deux impulsions 90° et 180°. L'amplitude A(TE ) du maximum du signal d'écho est telle que: _TE
  A(TE )
  = A(O)e
  (III. 1)
  T2
  où A(O) est l'amplitude du signal FID à l'instant initial t=O.
  66 Cette amplitude du maximum du signal d'écho, qui dépend du temps de relaxation transversale T2 de l'échantillon (et non de T;) , traduit le fait que l'influence des hétérogénéités de Bo est éliminée à l'instant TE' Plus de détails sur le mécanisme de la fonnation des signaux d'écho de spin se trouvent dans la littérature [Liang, 1999] et [Guinet, 1992] .
  Un exemple de signaux d'écho de spin obtenus avec notre système d'imagerie est illustré dans la figure III.3 qui donne la partie réelle de ce signal. Les excitations sont des impulsions sine 90° et 180°. Le temps d'écho TE étant de 20ms.
  : •• •• : • . •• 1•
  .
  :
  Impulsion 180°
  *1* FID initial
  Figure m.3
  Echo de spin
  Signal d'écho de spin obtenu avec une impulsion 90° suivie d'une impulsion 180°. Le temps d'écho TE est de 20ms.
  67 1/.2.2. Signal d'écho de gradient Un signal d'écho de gradient est produit par l'application adéquate d'un gradient. La chronologie de l'expérience conduisant à ce signal est donnée par la figure IlI.4.
  ____ ______________________________________ Impulsion RF
  ~n~
  t__).
  t
  Gradient
  )
  . A(0)  \
  ,A(O)e;;
  ****1",********** ,
  Signal
   '[ .
  t
  --------+-----~
  o
  )'
  T
  Figure m.4
  T
  T
  lllustration du signal d'écho de gradient
  Après l'application d'une impulsion RF, un gradient négatif de champ magnétique est appliqué suivant une direction donnée pendant une durée T. Ce même gradient est ensuit inversé et appliqué pendant une durée égale à 2T. On constate alors l'apparition d'un signal d'écho dont le maximum se situe au milieu du lobe positif du gradient. L'intervalle de temps ou temps d'écho, noté encore TE' qui sépare l'impulsion 90° du maximum du signal d'écho,
  68 est égale à la durée de l'application du gradient positif. L'amplitude A(TE ) du maximum du signal d'écho est telle que :
  _TE
  A(TE )
  = A(O)e
  (III.2)
  T;
  où A (0) est l'amplitude du signal FID à l'instant initial. L'amplitude du signal d'écho de gradient dépend du temps de relaxation effectif T2•• Le signal d'écho de gradient n'élimine donc pas, contrairement à celui d'écho de spin, l'influence des hétérogénéités du champ magnétique. Un signal d'écho de gradient acquis grâce à notre système d'IRM est présenté sur la figure III.5. Cet écho à été obtenu avec un une impulsion rectangulaire 90° et un temps d'écho
  TE égale à 12ms.
  Signal d'écho de gradient
  1
  Figure llI.5 Signal d'écho de gradient obtenu avec une impulsion rectangulaire 90° et un temps d'écho TE =12ms
  69 III.
  Construction de l'image en IRM Parmi les nombreuses techniques de reconstruction de l'image, nous allons présenter ici
  la technique d'imagerie par Spin-Warp qui est, de très loin, la plus utilisée, que ce soit pour l'imagerie 2D ou 3D [Guinet, 1992]. On se limite dans cet exposé au cas de l'imagerie 2D. Il s'agit alors de reconstruire l'image du plan d'une coupe dans l'échantillon à étudier. L'étape préliminaire de la réalisation d'une image consiste à sélectionner cette coupe. L'image de la densité protonique dans le plan de celle-ci peut être obtenue avec l'application adéquate de gradients dit de codage par la fréquence et de codage par la phase.
  111.1. Sélection d'une coupe dans l'échantillon: Cette sélection se fait à l'aide des impulsions dites sélectives. Elles sont produites par l'application simultanée d'une impulsion RF et de l'un des trois gradients du système d'imagerie. Ce gradient est appelé gradient de sélection. Nous allons, dans ce qui suit, expliciter la notion des impulsions sélectives. Nous choisissons le gradient selon l'axe Z, Gz , comme gradient de sélection et nous nous contentons d'exposer la procédure de sélection d'une coupe de cote z=o. Comme nous l'avons exposé dans le chapitre 1 (paragraphe IV) , lors de l'application d'un gradient selon l'axe Z, un plan, de cote Z, de l'échantillon se trouve soumis à un champ total B qui est la somme du champ magnétique Bo et du champ b = GzZ produit par le gradient Gz :
  (III. 3)
  Le champ magnétique total dans le plan de cote précession dans ce plan est donc
  z=o vaut donc
  Bo, et la fréquence de
  Jo = yBo . 27r
  La fréquence de précession, I(Z) , dans un plan de cote Z quelconque, est donc donnée par cette équation :
  (111.4)
  70 qui montre que cette fréquence varie linéairement avec la cote Z. (figure III.6).
  feZ)
  Z
  o Figure m.6
  Fréquence de précession en fonction de Z en présence d'un gradient de champ magnétique selon l'axe Z
  Cette variation de fréquence permet d'isoler (sélectionner) une certaine épaisseur dans l'échantillon. Pour cela, il faut appliquer, en plus du gradient, une impulsion d'excitation de fréquence centrale
  10
  et d'une certaine largeur fréquentielle l'!f. En effet, une telle impulsion
  n'agit efficacement que sur les moments magnétiques dont les fréquences de précession sont situées dans une bande fréquentielle de largueur l'!f centrée en fo' Compte tenu de la relation linéaire qui existe entre la fréquence de précession et Z, on en déduit que seuls les moments magnétiques situés dans une tranche de cote Z=O et d'épaisseur l'!Z (figure III.7), subissent l'action de l'impulsion. L'épaisseur de cette tranche (coupe) et la largueur fréquentielle l'!f de l'impulsion sont reliées par la relation :
  (III. 5)
  Cette relation montre que si Gz est égal à lmT/m, alors une impulsion de largueur l'!I égale à 426Hz sélectionne une coupe de lem d'épaisseur. La même impulsion, si Gz vaut lOmT/m, permet de sélectionner une coupe de lmm d'épaisseur.
  71 feZ)
  :<
  ): !U,
  Figure ill.7
  Sélection d'une coupe d'épaisseur !U, et de cote Z=O, par une impulsion de fréquence 10 et de largeur fréquentielle /).f .
  Le signal de résonance reçu après une impulsion sélective provient uniquement de la coupe sélectionnée. Il ne contient que les informations relatives à cette coupe et constitue une donnée globale sur celle-ci. Il s'agit maintenant de construire, à partir de ce signal global, l'image de la distribution des protons dans le plan de la coupe (plan XOY). Ceci peut être accompli avec l'utilisation des gradients Gx et Gy.
  111.2. Codage de fréquence et codage par la phase : reconstruction de l'image Après qu'un signal a été activé par une impulsion sélective, les informations spatiales, relatives à la distribution de la densité de protons dans le plan de la coupe, peuvent être codées dans ce signal pendant la période de précession libre. Compte tenu du fait que le signal RMN peut être représenté par une exponentielle complexe, nous avons essentiellement deux possibilités pour coder des informations spatiales dans ce signal: codage par la fréquence et codage par la phase du signal [Liang, 1999]. 111.2.1 Codage par la fréquence
  Pendant la période de l'acquisition du signal, l'application un gradient de champ magnétique, Gx , selon l'axe X, entraîne une dispersion linéaire des fréquences de résonance dans la direction OX Ainsi, la fréquence de précession à l'abscisse X s'écrit :
  72 f(X)
  = fa + ,.cx x
  (III. 6)
  2n
  Le signal élémentaire recueilli provenant d'un élément dX d'abscisse X, a pour expression [Liang, 1999] :
  dS(X,t)
  = p(X)e i2Trf (fo +Gx X)t dX
  (III.7)
  où p(X) est la densité de protons de l'élément dX Ce signal est dit codé par la fréquence car sa fréquence dépend linéairement de sa localisation spatiale suivant l'axe X. Le gradient Gx s' appelle ainsi gradient de codage par la fréquence. Le signal total, qUI est la somme de tous les SIgnaux élémentaires, a donc pour expreSSIOn :
  S(t)
  J p(X )e
  =
  i 2Trf
  (fo+ Gx
  X)t dX
  (III.S)
  échantillon
  Après démodulation, i.e. translation de la fréquence fa, et en posant K x
  = ,.cxt , cette
  expression est réécrite comme suit:
  S(K x)
  =
  J p(X)e
  i2
  x XdX
  (III.9)
  XdKX
  (III. 10)
  Jd(
  échantillon
  d'où: p(X) =
  J S(K)e-
  i2 Jd(
  x
  échantillon
  La transformée de Fourier inverse de S(K x ) permet donc d'avoir la distribution de la densité de protons le long de l'axe X. Il est important de noter que le gradient de codage par la fréquence est appliqué pendant la durée de l'acquisition du signal S(t). Nous choisissons désormais ce gradient de codage par la fréquence pour coder l'espace selon l'axe X
  73 111.2.2 Codage par la phase
  Cette méthode consiste à manipuler la phase initiale du signal RMN acquis. Pour comprendre la notion du codage par la phase, imaginons qu'on applique un gradient,
  Gy,
  selon OYpendant une durée Ty avant l'acquisition du signal RMN. Le signal provenant d'un élément dY, sous l'influence de ce gradient, sera donné par cette équation :
  p(Y)e i2 1l)'(fo+GYY )1 dY dS(Y,t)
  =
  (IIL11 )
  { p(Y)e i2 1l)'(G yy)Tyei2trtfol dY
  Cette relation montre que, pendant l'intervalle de temps 0:$ t:$ Ty , le signal local est codé par la fréquence. Le résultat de ce codage par la fréquence est que les signaux de différentes positions Y vont acquérir des phases différentes après un temps Ty • Si on utilise cet intervalle de temps, 0 :$ t :$ Ty , comme une période de p réparation à l'acquisition du signal, celui-ci va donc avoir une phase initiale donnée par la relation [Liang, 1999]
  <D(Y ) = t<IyYTy
  qui peut se réécrire, en posant Ky
  (IlL 12)
  = t<IyTy , comme suit: (IlL 13)
  111.2.3 Association du codage par la fréquence et du codage par la phase
  Pendant la période de préparation Ty , le gradient de codage par la phase est appliqué. A l'issue de cette période, i.e. à partir de l'instant t;;:: Tf' on applique le gradient G x signal élémentaire recueilli d'un élément de volume dXdY sera donnée par :
  dS(X,Y,t)
  = p(X,Y)ei2nrGYYTy e i2 1l)'(fo+GxX)1 dXdY
  où p(X,Y) est la densité de protons du volume élémentaire dXdY.
  Après démodulation, l'expression III. 14 devient : 74
  (ilL 14)
  .
  Le dS (X, Y, t)
  = p( X, Y)e i2 ll)'G yYTyei2TrY {GxX)t dXdY
  (IILl5)
  ce qui conduit, en tenant compte des expressions de K x et Ky, à :
  (III.l6)
  Le signal total recueilli pendant l'application du gradient de codage par la fréquence est donc donné par cette équation :
  S(Kx,Ky)
  =
  ffp(X , Y)ei2 1r{KxX+Kyy)dXdY
  (IILl7)
  échantillon
  Pour déterminer l'expression de p(X,Y) et donc l' image du plan de la coupe, il suffit d'effectuer une transformé de Fourier en deux dimensions (2D) inverse, soit:
  p(X,Y ) =
  ffS( K x ,K y)e-i21r{KxX+Kyy )dK XdKy
  (IIU 8)
  échantillon
  Dans cette équation, qui constitue l'équation de base de l'IRM par transformée de Fourier, l'ensemble de points (Kx,Ky) définit un plan appelé plan de Fourier. C'est, en effet, l'espace réciproque du plan de l'image (le plan réel) défini par l'ensembles des points (X, Y) [Liang, 1999]. La tâche qui reste donc à accomplir pour construire une Image est de collecter suffisamment de données pour remplir l'espace de Fourier. Pour cela, une stratégie conventionnelle consiste à remplir ligne par ligne} cet espace de Fourier. C'est la méthode de base de l'acquisition d'une image.
  1II.2.4.Acquisition de l'image, balayage du plan de Fourier:
  1 Le balayage ligne par ligne ne constitue que le principe de base. D 'autres trajectoires du balayage du plan de Fourier peuvent être réalisées. Tout dépend de la séquence d'imagerie choisie
  75 111.2.4.1. Principe du balayage
  En pratique, on fait l'acquisition des points discrets dans le plan de Fourier. Supposons, par exemple qu'on cherche à acquérir une image (une matrice de points) de N y lignes et de N x colonnes. Une ligne de cette matrice correspond à une valeur donnée de Ky
  = yGyTy
  (donc des
  valeurs données de Gy et de Ty ). Avant l'acquisition du signal, on applique un gradient Gy de codage par la phase pendant une durée Ty • Une ligne de la matrice a été donc fixée. Les valeurs K x
  = yG xt ,
  qui sont en nombre N x , de cette ligne sont acquises en prélevant N x
  échantillons du signal pendant la durée de l'acquisition où le gradient de codage par la fréquence Gx
  est appliqué. L'échantillonnage du signal se fait avec une période
  d'échantillonnage
  I:. La durée d'acquisition du signal est donc T acq = N XTe'
  Pour remplir une deuxième ligne, il faut modifier la valeur Ky
  = yGyTy
  avant
  l'acquisition d'un nouveau signal avec le même gradient de codage par la fréquence G x . La modification de la valeur Ky peut être effectué en modifiant soit Ty, soit Gy. Dans tout ce qui suit, nous considérons que c'est la valeur de Gy qui sera modifiée, d'une ligne à l'autre, avec un pas fixe qu'on notera désormais oGy. Ainsi, pour acquérir les N y lignes, il faudra répéter N y fois la même procédure d'acquisition du signal, avec à chaque fois une valeur différente du gradient de codage par la phase. Par exemple, l'acquisition de la ligne numéro m se fait en appliquant un gradient de phase moG y (figure 111.8). Dans cette figure, les valeurs K x de la colonne numéro n correspondent l'échantillon du signal pris à l'instant nTe avec des valeurs différentes de Ky. La transformée de Fourier 2D discrète de la matrice (N x x N y) de données ainsi acquise donne alors l'image de la densité de protons recherchée.
  76 ~t.
  KY
  Ligne « m >.:> Ky =rm8G yT} Kx
  
  0
  t Colonne « n »
  Kx
  Figure m.8
  = rGx(nTe )
  Les valeurs de K x et de Ky correspondant à une ligne m et à une colonne n de l'espace de Fourier.
  111.2.4.2. Considérations pratiques, paramètres d'acquisition
  Le nombre N x d'échantillons prélevés au cours de la période d'acquisition du signal détermine le nombre de colonnes de largueur !lX découpant l'espace réel (figure III.9). Ces N x colonnes occupent une région de l'espace réel limitée dans la direction X La dimension Lx de cette région est:
  (III. 19)
  77 ,--.,--.,-,---,,,,,,"","-.,--.,-,-,-,---, 
  x
  )
  (
  Figureill.9
  Les voxels dans l'espace réel.
  De même, le nombre N y détermine le nombre de lignes de largueur
  ~Y,
  donc le
  nombre de voxels de chaque colonne. Ces N y voxels délimitent, dans la direction Y, une région de longueur Ly de l'espace réel:
  (lII.20)
  Par conséquent, Lx et Ly déterminent une région dans l'espace réel qu'on appelle champ d'imagerie (Field ofView (FOV)) [Guinet, 1992]. En général, Pour une matrice de taille donnée (N x x N y), l'utilisateur peut choisir ce champ d'imagerie. La résolution par pixel de l'image sera donc: M
  = Lx (ou ~y = Ly ). Nx
  Ny
  Cependant, d'autres considérations primordiales doivent être prises en compte. Cellesci comprennent, d'abord, la fréquence d'échantillonnage, Fe' du signal lors de l'acquisition qui doit, bien évidemment, respecter le théorème de l'échantillonnage de Shannon.
  78 La fréquence la plus élevée dans le signal est:
  r
  (III.21)
  fmax. =-GXXmax. 27r , X ou max.
  Lx =2
  La fréquence d'échantillonnage2 doit donc vérifier la relation suivante:
  (III.22) ou encore:
  (III.23)
  De façon identique, la dimension
  ~Y
  des voxels suivant l'axe Y, le pas de gradient de
  phase oGy, la durée d'application de ce gradient Ty, et le nombre N y doivent vérifier la relation d'échantillonnage spatiale selon Ky [Guinet, 1992] :
  (III. 24)
  A partir des deux relations III.23 et III.24, nous devons prendre en compte, dans notre système d'IRM, le fait que Lx et Ly ne peuvent pas être abaissés en deça d'un certain seuit3 ., pour les deux raisons suivantes: •
  A Fe fixé, puisque, pour un appareil donné d'IRM, le gradient Gx ne peut pas dépasser une certaine valeur maximale caractéristique de l'appareil, le choix de Lx n'est accepté que si la valeur de Lx est supérieure à une valeur minimale, inversement proportionnelle à la valeur maximale de Gx (relation III.23).
  2
  Pour notre récepteur numérique la fréquence
  Fe
  présente la fréquence d'échantillonnage finale dans la bande
  de base. Nous avons choisi de la noter ainsi pour la distinguer de du symbole l e utilisé dans le chapitre II. Il va aussi de soi que le champ d'imagerie ne peut pas être abaissé en deça des dimensions de l'échantillon à imager car cela produit des artéfacts de repliement [Guinet 1992].
  3
  79 •
  Pour une durée Ty fixée, la valeur de Ly choisie par l'utilisateur permet à partir de l'équation III.24 de déterminer la valeur de oGy et, donc, toutes les valeurs successives du gradient selon l'axe Y pour les N y acquisitions. Aucune de ces valeurs ne devant dépasser la valeur maximale du gradient applicable selon l'axe Y. Le choix de Ly n'est accepté que si la valeur de Ly est supérieure à une valeur minimale, liée à la valeur maximale du gradient selon Y.
  Nous avons effectué sur notre système, les mesures nécessaires pour déterminer les limites des gradients que peuvent créer les bobines de gradients de l'électroaimant. Grâce à ces mesures, nous avons pu établir que les valeurs maximales des gradients selon les trois axes sont identiques :
  Dans le paragraphe suivant, nous allons expliquer l'implantation de la séquence d'écho de gradient sur notre système.
  111.3. Les premières images obtenues avec le système numérique d'IRM par la séquence d'écho de gradient: 111.3.1. Réglage de l'angle de basculement de l'impulsion
  Cette procédure vise à déterminer l'amplitude (donc la puissance) et la durée de l'impulsion permettant d'obtenir l'angle de basculement désiré pour la séquence. Ce réglage doit être effectué pour chaque échantillon avant de procéder à l'acquisition de l'image. Nous avons vu, au chapitre l, que l'angle de basculement dépend du produit de l'amplitude de l'impulsion par sa durée. Dans une séquence d'imagerie, pour un gradient de sélection donné, la durée
  l',
  qui est proportionnelle à l'inverse de la largeur fréquentielle de
  l'impulsion, est fixée par l'épaisseur de la coupe qu'on souhaite obtenir. L'angle de basculement doit donc forcément être réglé en agissant sur l'amplitude de l'impulsion. Pour ce faire, il suffit d'augmenter progressivement la puissance de cette dernière jusqu'à l'obtention du signal RMN maximal. L'angle de basculement correspondant est donc de 90°. Dans notre système d'imagerie, une séquence d'écho de gradient de temps de répétition égal à 1 seconde est utilisé pour ajuster l'angle de basculement. La procédure de
  80 l'ajustement est la suivante4
  :
  l'échantillon à imager est placé dans l'antenne et le signal
  d'écho de gradient est observé à l'oscilloscope. La puissance de l'impulsion est modifiée jusqu'à ce qu'on observe un maximum du signal d'écho. L'angle de basculement de l'impulsion appliqué est donc de 90°. De même, pour obtenir une impulsion 180°, il suffit d'augmenter la puissance de l'impulsion pour ne détecter aucun signal d'écho. Pour calculer la puissance nécessaire à obtenir un angle de basculement quelconque
  a, une solution consiste à procéder par un simple règle de trois à partir de l'énergie nécessaire pour créer une impulsion 90°. 111.3.2. La séquence d'écho de gradient
  Le chrono gramme d'un cycle de la séquence d'écho de gradient, que nous avons implanté sur notre système, est représenté dans la figure Ill. 1O. L'impulsion RF est appliquée simultanément avec le gradient de sélection de la coupe
  Gz . Le gradient de codage par la phase, Gy , est ensuite appliqué pendant une durée Ty. Le lobe négatif du gradient de codage par la. fréquence G x est appliqué pendant la durée de l'application du Gy. L'acquisition du signal est effectuée pendant la durée Tacq ' Le signal présente un maximum à l'instant TE (temps d'écho). Le cycle de la séquence est répété N y fois avec un temps de répétition TR • Seule la valeur du gradient de phase est modifiée d'une répétition à l' autre.
  Dans le cas d'une impulsion sélective, lors de l'ajustement de l'angle de basculement, il est évidemment nécessaire d'appliquer le gradient de sélection de coupe.
  4
  81 Impulsion
  ~------------~m~ i ,
  ,
  ,
  1
  Gi1
  ! '
  Gy
  1
  1
  <
  )
  Ty Gx
  L Signal
  Acquisition Figure m.tO
  La séquence d'écho de gradient.
  La séquence a été programmée en langage assembleur de la carte DSP. Celle-ci contrôle tous les événements temporels de la séquence. Les signaux acquis sont stockés dans la mémoire du DSP. Les premières images, d'une coupe d'épaisseur de lcm obtenue par une impulsion sélective rectangulaire 90°, ont été réalisées avec une fréquence d'échantillonnage et une durée d'application du gradient de phase fixées à Fe = 10kHz et Ty = 3ms . Le temps d'écho
  82 et les temps de répétition sont TE = 6ms TR = 1.2s. Nous avons choisi une matrice d'acquisition de N X x N y= 64x64 et un champ d'imagerie de LXxLy = 12x12cm. La résolution obtenue dans ces conditions est d'environ 1.8mm/pixel. Un exemple des images ainsi obtenues, sur un tube cylindrique d'eau déminéralisée, d'environ 3cm de diamètre, est illustré par la figure IIL11.
  10
  20
  30
  .0
  '"
  60
  Figure m.u Images 64x64 pixels obtenues avec la séquence d'écho de gradient. A gauche, l'image d'une coupe dans un échantillon cylindrique d'eau déminéralisée de 3cm de diamètre. A droite le même tube quand on insère au centre objet cylindrique de 1cm de diamètre. L'excitation est une impulsion rectangulaire 90°. L'épaisseur de la coupe est 1cm. Le temps de répétition TR = 1.2s , et le
  temps d'écho TE = 6ms. La résolution étant de 1.8mm1pixel. Ces images ont été obtenues avec 8 accumulations et un temps total d'acquisition de 10 minutes.
  IV. •
  Quelques travaux en cours La méthode de la bande passante utilisée pour la mesure de la fréquence de résonance, du facteur de qualité et pour 1'accord en impédance de l'antenne présente les inconvénients de manquer de précision et de la difficulté de raccorder la fréquence de résonance et l'impédance de l'antenne pour chaque échantillon dans les conditions d'acquisition. Nous travaillons actuellement sur le développement d'un module permettant, par l'utilisation d'un coupleur de puissance RF, d'accorder l'antenne. Les modifications de la courbe de résonance, apportées par la boucle d'accord de l'antenne, seront affichée en temps réel sur l'écran de l'ordinateur ce qui permet un réglage très fin de l'accord pour chaque échantillon.
  •
  Implantation de la séquence d'écho de spin.
  83 •
  Interface homme-machine: nous travaillons également sur le développement d'une interface graphique conviviale pour notre système d'imagerie. Cette interface utilise un environnement Windows NT. Elle permettra à l'utilisateur de définir facilement tous les paramètres d'acquisition de l'image, en générant un fichier pouvant être interprété par le DSP; de visualiser et de sauvegarder toutes les informations relatives à l'image. Cette interface permettra aussi au développeur de séquences d'IRM d'implanter de nouvelles séquences d'imagerie sur notre systèmes.
  •
  Après une étape de validation et de comparaison des images acquises avec celles obtenues avec un imageur utilisant une électronique analogique, nous envisageons l'implantation des séquences ultra-rapides pour l'imagerie en temps réel. Ceci peut être autorisé sur notre système grâce au bon RSB. En effet, le RSB est un des facteurs déterminant du temps d'acquisition d'une image. Ce temps est souvent long (quelques minutes) à bas champ à cause des opérations de moyennage du signal qui sont nécessaires pour améliorer le RSB. S'il est élevé, ce dernier permet donc de limiter le nombre de moyennages et par conséquent de réduire le temps d' imagerie.
  En conclusion à ce chapitre, nous avons pu valider le développement numérique d'une IRM bas champ. Les premières images esquissées montre la faisabilité de notre démarche de numérisation. Outre le coût peu élevé, par rapport à un système analogique moins performant en terme de RSB et de flexibilité, la solution numérique adoptée permet de miniaturiser les systèmes d'IRM. En effet, notre système comporte, hors aimant, un PC où sont enfichées trois cartes électroniques (la carte DSP superposée avec la carte récepteur numérique et avec celle du synthétiseur numérique) et un amplificateur RF de puissance peu volumineux. Ajoutons finalement que, grâce au récepteur et au synthétiseur numériques, l'électronique mise en place et destinée
  ICI
  au bas champ, peut être, à l'exception de
  l'amplificateur RF de puissance, intégrée pour fonctionner dans des systèmes d'imagerie à hauts champs magnétiques pouvant atteindre 3.5T.
  84 Chapitre IV
  Identification des systèmes RMN par filtrage adaptatif non linéaire : Modèle de série de Volterra
  85 1.
  Introduction La plupart des concepts de la spectroscopie et de l'imagerie par résonance magnétique
  nucléaire sont basés sur la considération de systèmes linéaires ou approximativement linéaires pour lesquels un traitement mathématique simple et élégant est possible, et pour lesquels nous avons en quelque sorte une compréhension intuitive. Néanmoins, un système RMN (échantillon à étudier) est fondamentalement non linéaire, et l'application de la théorie de la réponse linéaire nécessite quelques précautions [Ernst, 1989] . Pour des raisons de meilleure représentation des systèmes RMN, on sera donc conduit à envisager la prise en compte des caractéristiques non linéaires de tels systèmes.
  Contrairement à un système linéaire qui est complètement identifié par sa réponse impulsionnelle, l'élaboration d'une méthode générale d'identification pour les systèmes non linéaires est loin d' être triviale. En effet, la méthode à adopter pour l'identification non linéaire dépend fortement du problème à traiter et de la nature des applications ultérieures tout en prenant en considération que l'amélioration des performances du modèle linéaire doit justifier la complexité du calcul accompagnant généralement les problèmes non linéaires. D' une façon générale, les propriétés d'un système, considéré comme une « boîte noire », sont formellement représentées par un opérateur cI> liant l'entrée et la sortie. Pour le système RMN, l'entrée x(t) correspond à l'excitation RF, ce qui conduit en sortie au signal FID observable d(t) (figure. 1).
  x(t)
  Figure IV.l
  
  d(t)
  Boîte noire (système RMN) cI>
  
  La relation entrée-sortie d'un système donné est représentée par un opérateur cI>.
  L'objectif de l'identification consiste à établir un développement de l'opérateur cI> du système. Après une brève analyse de la non linéarité associée aux systèmes RMN, et dans la perspective d'optimiser le traitement et l'analyse du signal RMN, nous allons introduire dans ce chapitre un modèle non linéaire basé sur la description de la relation entrée-sortie du
  86 système RMN par une série de Volterra. Les trois premiers noyaux de celle-ci seront identifiés à l'aide de techniques de filtrage adaptatif non linéaire. Le modèle proposé sera validé sur des signaux RMN obtenus grâce au simulateur décrit dans le premier chapitre.
  Il.
  Non linéarité de systèmes RMN et nécessité de l'identification non linéaire Une expérience RMN consiste donc à appliquer une excitation BI (t) à l'entrée d'une
  «boîte noire» contenant un échantillon afin d'observer la réponse (FID) à cette excitation. L'objectif étant bien sûr d'étudier certaines propriétés de cet échantillon. La réponse n'est généralement observable qu'après la fin de l'excitation, c'est-à-dire à partir de l'instant où BI (t) devient nul. La relation entre l'aimantation, M xy (t) pendant la durée de la précession
  libre et l'excitation s'exprime, dans le repère tournant, par l'équation suivante:
  M Xy (t) =Mo sina e Yr2
  e-il!.(i)of
  =Mosin( f rBt(t')dt'
  ) e
  Yr2
  e-il!.(i)of
  (IV. 1)
  o
  Cette relation confirme le fait que la réponse d'un système RMN est une fonction non linéaire de l'impulsion d'excitation. TI suffit pour s'en convaincre de constater que si une impulsion
  f agissant sur un système RMN produit un certain signal d(t), une impulsion 7t ne
  produit pas un signal 2d(t), comme le ferait un système linéaire, mais un signal nul. Cependant, si l'impulsion est de durée et/ou d'amplitude suffisamment faible pour ne perturber que légèrement l'échantillon, il est possible de montrer que le système se comporte comme un système linéaire. Cette approximation linéaire est valable tant que l'angle de basculement reste inférieur à 30° [Hoult, 1979]. Pour des angles de basculement plus élevés, le système est conduit dans la région non linéaire où le traitement linéaire d'une réponse intrinsèquement non linéaire peut induire des erreurs aussi bien de traitement que d'analyse ultérieure du signal (distorsion du spectre par exemple) [Blümich, 1982]. Des approches d'identification non linéaires s'avéreront donc indispensables pour assurer le traitement correct des signaux. Pour ce faire, nous avons intérêt à exprimer la relation (IV.1) sous une forme adaptée aux techniques de traitement numérique.
  87 III.
  La série de Volterra appliquée à l'identification de systèmes RMN Introduite au début du 20ème siècle, la série de Volterra peut théoriquement décrire une
  classe très large de systèmes non linéaires. Elle se présente, généralement, sous la forme d'une somme infinie de produits de convolution multidimensionnels [Kaulisch, 1996], [Asfour-JMR, 2000] :
  =
  d(t)
  f p=O
  r r . .r
  hlrI' r 2 ,
  ••• r
  p) x(t-r l ) x(t-r 2 )***x(t - r p) dr l dr 2 ***dr p
  (IV.2)
  Les fonctions temporelles multidimensionnelles h p(ri' r 2'  , r p)' souvent appelées noyaux de Volterra, caractérisent la non-linéarité. Cependant, en l'absence d'une méthode générale pour l'identification de ces noyaux, cette série n'a, en général, qu'un intérêt pratique très limité, et malgré sa longue histoire et sa popularité dans les études théoriques, peu de recherches ont tenté d'identifier des systèmes pratiques en se basant sur cette représentation. Ceci est attribuable aux difficultés associées à l'identification des noyaux de Volterra [Fakouri, 1980], [Lawrence, 1981], [Koh, 1985] et [Biglieri, 1982].
  111.1. Forme discrète de la série de Volterra: Nous développons dans ce qui suit une méthode d'identification basée sur un filtrage adaptatif non linéaire. Cette méthode considère la forme discrète de la série de Volterra qui s'exprime alors par la relation suivante [Mathews, 1991] :
  co
  d(n)
  =
  'I hlm)x(n - ml) +
  co
  co
  co
  n~=O
  m2=O
  mp =O
  'I 'I 
  co
  +
  co
  'I 'I hlml' m
  'Ihlml' m 2 ,
  ••• ,
  2)
  x(n - ml )x(n - ml) + 
  mp )x(n-m l )x(n-m 2 )
  88
  •••
  x(n-m p )+
  (IV.3) où x(n) et d(n) sont respectivement l'entrée et la sortie à l'instant n du système discret causal et non linéaire. En tronquant cette série à un nombre fini de termes, on obtient une approximation à l'ordre p de la forme [Asfour, 2000]:
  N-'
  d(n)
  =
  N-' N-'
  Ihlm,)x(n - m,) +
  +
  N-' N-'
  I I  n/I = O
  n/2 =O
  l
  Ihim" m 2 )x(n - m,)x(n - m) +  (IV.4)
  N-'
  Ih/m" m 2 ,
  •••,
  mp) x(n-m,)x(n-m 2 )
  x(n-m p)
  •• •
  11I p =O
  où N est la mémoire du système et hp (m" m2 , •• •, mp) est le noyau de Volterra d' ordre p .
  L'identification du système se résumera alors au calcul des coefficients des p noyaux de Volterra qui décrivent à leur tour le comportement du système. La précision dépend de l'ordre de l'expansion p et de la mémoire N choisis.
  111.2. Identification par série de Volterra de second ordre
  111.2.1. Méthode et schéma général d'identification
  Nous considérons le développement à l'ordre 2 de la série de Volterra (p=2), soit donc:
  d(n)
  =
  N-'
  l
  N-'
  hlm,) x(n - m,) +
  N-'
  l l
  him" m 2 ) x(n-m, )x(n-m 2 )
  (IV.S)
  Le diagramme général correspondant à l'identification d'un système RMN est représenté sur la figure IV.2
  89 SystèmeRMN r
  Signal de référence
  
  x( n)
  den)
  Filtre de Volterra signa1 d'entrée
  ~
  I~
  1
  retard
  Filtre linéaire
  -+ .i
  r--"-
  È~
  / Filtre quadratique
  -+
  d(n) \,-
  .    ~ :
  ~
  r--
  erreur ~
  e(n) = den) - den) ~
  N-l
  ~
  d(n) = L~(ml : n)x(n-m1) +
  Figure IV.2
  N-l N-l
  I
  ~
  I~(ml,m2 :n)x(n-m1)x(n-m 2)
  Schéma général d'identification par filtre adaptatif de Volterra de 2ème ordre.
  Le filtre adaptatif de Volterra est donc constitué de l'association parallèle d'un filtre linéaire et d'un filtre quadratique. Le signal d'entrée x(n), i.e. l'impulsion d'excitation, est appliqué au système RMN à étudier dont la réponse d(n) est désignée dans le contexte de l'identification par le terme signal de référence. Les paramètres du filtre adaptatif doivent être choisis de manière à ce que le modèle obtenu soit le plus proche possible du système à analyser. Le filtre de Volterra va donc essayer de calculer, une estimation
  den) de la réponse
  désirée en utilisant une série de Volterra tronquée à l'ordre 2. Cette estimation est donnée à l'instant n par :
  N-I
  d(n)
  =
  L hlm l.' n) x(n - ml)
  N-l
  +
  A
  N-l
  L L hlmp m
  2 .'
  n) x(n-m l )x(n-m 2 )
  (IV.6)
  A
  Dans cette équation, hl (ml : n) et h2 (ml ,m 2 : n) , avec 0 ~ ml , m 2
  ~
  N -1, sont
  respectivement les coefficients du filtre linéaire et du filtre quadratique à l'instant n. Le
  90 problème posé consiste à déterminer, à chaque n, l'ensemble de ces coefficients qui minimise une fonction coût du signal d'erreur e(n) entre d(n) et den) donné par:
  e( n)
  = d (n ) -
  A
  (IV.7)
  d (n )
  En pratique, le calcul des coefficients peut être mis en oeuvre d'une manière récursive. Des algorithmes récursifs ont pour objectif de calculer les coefficients à l'instant n+ 1 à partir de ceux à l'instant n et ce avec une initialisation correcte des valeurs de départ (n=O). Notons qu'en résonance magnétique, l'entrée et la sortie «observable» du système ne sont jamais simultanées (la sortie n'étant pas observable pendant l'excitation). Ceci doit être pris en compte dans le modèle proposé. Pour ce faire, nous introduisons un retard temporel L1 sur le signal d'entrée du filtre (figure IV.2). La valeur de ce retard dépend essentiellement de la mémoire du système et, bien entendu, du modèle à estimer [Widrow, 1975] [Raoof, 1994]. Nous considérons également que le signal RMN est nul pendant l'excitation.
  111.2.2. Algorithme d'adaptation Nous avons adopté l'algorithme de gradient classique pour résoudre le problème d'identification des systèmes RMN. Ce choix est justifié par la simplicité de cet algorithme et la souplesse de sa mise en oeuvre. L'algorithme est basé sur la minimisation de l'erreur quadratique e 2 (n) à chaque instant n. Les coefficients hl (m 1 : n) et h2 (ml' m2 : n) , avec 0 ~ ml' m 2
  ::;;
  N -1 , du filtre sont
  mis à jour par les équations suivantes [Mathews, 1991]:
  (IV. 8)
  (IV.9)
  91 Les scalaires 0 1 et O2 , appelés pas d'adaptation, contrôlent la stabilité de l'algorithme adaptatif et sa rapidité de convergence. En pratique, la mise en oeuvre de cet algorithme s'effectue en adoptant des notations vectorielles permettant, hormis des calculs plus rapides, une meilleure analyse des performances. A l'aide de ces notations, il est possible d'écrire l'estimation
  den)
  sous la
  forme suivante :
  (IV. 10)
  où HI(n) est le vecteur des coefficients du filtre linéaire, et H 2(n) est le vecteur des coefficients du filtre quadratique à l'instant n :
  Hln)
  =
  [hlO: n), hl1: n),  , hl(N - 1 n) ] T
  Hin)
  =
  [hiO,O: n), hi0,1: n),  ., hlO, N -1: n), hil,l: n), hi 1,2: n),  , h 2(N-1, N-1: n)] T
  (IV. 11)
  (IV. 12)
  X I(n) et X 2 (n) sont respectivement les vecteurs d'entrée du filtre linéaire et du filtre
  quadratique. Ils sont définis à l'instant n par:
  Xln)
  = [x(n), x(n -1),  , x(n -
  Xin)
  = lxYn), x(n)x(n -1),
  N + 1)] T
  (IV. 13)
   ,x(n)x(n - N + 1), x 2(n -1), x(n -l)x(n - 2),  , x 2(n - N + l)J T (IV. 14)
  Les équations de mise à jour s'écrivent donc .:
  e(n)
  = den) - Hr (n)X l (n) - HJ (n)X 2(n)
  (IV.1S)
  H l (n + 1) = Hl (n) + OIXl (n)e(n)
  (IV. 16)
  H 2(n + 1) = H 2(n) + 02X2 (n)e(n)
  (IV. 17)
  92 L'application de cet algorithme à l'identification des système RMN nécessite, bien entendu, un choix adéquat de ses paramètres caractéristiques à savoir les pas d'adaptation 0 1 , 8 2 et la mémoire du système N C 'est ce que nous analyserons par la suite.
  111. 2.2.1.
  Choix des pas d'adaptation pour les signaux RMN
  Dans le cas général d'un filtre quadratique de Volterra, la stabilité de l' algorithme est assurée si les pas d ' adaptation 8 1 et 8 2 sont choisis tels que [Bellanger, 1989] :
  (IV.18)
  (IV.19)
  où O'~ est la puissance du signal d'entrée x(n). Ces deux relations laissent apparaître que la stabilité de l'algorithme adaptatif dépend fortement de la puissance du signal d'entrée En résonance magnétique, les signaux d'entrée et de sortie du système RMN sont non stationnaires. Par conséquent, l'utilisation de valeurs constantes pour les pas d'adaptation peut conduire à des instabilités compte tenu de la nature impulsive des signaux. Pour remédier à ce problème, nous avons introduit des pas d'adaptation variables de la forme:
  bl(n)
  bln)
  où
  PI
  et
  /32
  = =
  /31
  (IV.20)
  N(J~(n)
  /32
  (IV.21)
  (NO'~ (n)) 2
  sont des facteurs positifs et O'~ (n) est une estimation de la puissance du signal
  d'entrée à l'instant n. L'estimation la plus directe de O'~ (n) est donnée par :
  (IV.22)
  93 avec (J'~ une constante positive pennettant d'éviter une division par zéro. La fenêtre d'observation temporelle No définit la durée sur laquelle le signal peut être considéré comme stationnaire.
  111.2.2.2.
  Choix de la mémoire du système N
  Les coefficients des deux noyaux de Volterra calculés * à l'aide de l'algorithme adaptatif sont la meilleure estimation du modèle du système RMN au sens de l'erreur quadratique moyenne (EQM). La précision de cette estimation est contrôlée par la mémoire du système N comme le montre la figure IV.5. Dans cette figure, sont reportées les variations de l'erreur quadratique moyenne en fonction de N. Cette erreur est calculée en utilisant la différence quadratique entre le signal de référence et le signal à la sortie du filtre en réponse au même signal d'excitation x(n). Le signal à la sortie du filtre de Volterra a été reconstruit en appliquant l'équation IV.5 et ce en utilisant les coefficients issus de l'identification. Nous pouvons voir que l'erreur quadratique moyenne est une fonction décroissante de N et atteint une valeur asymptotique pour les valeurs de N supérieures à 64. Cette valeur
  constitue donc un choix optimal pour la mémoire du système car elle pennet d'assurer un compromis entre le nombre de coefficients à estimer et la précision de l'estimateur [AsfourJMR,2002].
  94 10
  X 10-6
  EQM
  9
  8 7
  6 5
  4 3
  ~
  2 1 0
  20
  40
  60
  *
  *
  80
  100
  * 120
  140
  N
  Figure IV.3 Variations de l'erreur quadratique moyenne (EQM) entre le signal de référence et la sortie du filtre de Volterra en fonction de la mémoire du système N.
  111.2.3. Résultats de l'identification et discussions 111.2.3.1.
  Résultats
  Nous avons implanté l'algorithme d'adaptation sous MATLAB et nous l'avons employé dans un premier temps pour le choix de la mémoire du système comme nous l'avons vu ci-dessus. Pour illustrer les résultats de l'identification, nous utilisons des signaux RMN obtenus par la simulation. Nous précisons, préalablement, que le nombre de coefficients à estimer est de N 2 + N. Pour une mémoire système N=64, celui-ci vaut donc 4160. En prenant en compte le fait que la matrice donnant le deuxième noyau de Volterra est symétrique [Bellanger, 1989], il est donc suffisant d'effectuer l'estimation de la moitié de ses coefficients. Pour permettre cette estimation, il faut disposer d'un nombre suffisant d'échantillons du signal de référence. Ceci est réalisé par l'algorithme d'adaptation qui utilise les signaux
  95 correspondant à la répétition de la même impulsion. Cette même procédure sera décrite en détail dans le cinquième chapitre dans lequel nous traitons de signaux RMN réels. Dans le cadre du présent chapitre, les signaux utilisés par l'algorithme sont, bien évidemment, parfaitement identiques car il s'agit de signaux simulés. En se basant sur cette procédure, nous donnons ici les résultats de l'identification de l'objet ID décrit au chapitre I. Le signal d'entrée x(n) étant une impulsion sine 90° de lms de durée et le signal de référence den) étant présenté par la figure (1.8). L'estimation des deux noyaux de Volterra a été effectuée en fixant la mémoire du système à N=64. Dans ce cas, les valeurs des paramètres permettant de minimiser l'erreur quadratique moyenne (EQM) entre le signal de référence et la sortie du filtre de Volterra sont:
  f3l
  = 0.25,
  132
  =0.0025, CT; =0.08 , No = 20, et le retard Ll = 200 .
  Dans la figure IV.4 sont superposés le signal de référence et le signal à la sortie du filtre de Volterra.
  0.8
  0.6
  0.4
  02
  400
  450
  500
  550
  600
  650
  700
  n Figure IV.4 Le signal de référence (traits continus) et le signal de sortie du filtre de Volterra (en pointillés) pour une mémoire de système de 64. Les premiers 400 échantillons correspondant à la durée de l'impulsion ne sont pas présentés.
  96 Les coefficients du filtre linéaire (noyau de Volterra de premier ordre) et du filtre quadratique (noyau de Volterra de second ordre) correspondant sont illustrés par la figure IV.5. Dans cette figure, nous pouvons voir que le premier noyau de Volterra est une fonction s'apparentant au signal FID [Ernst, 1989]. La partie réelle est la partie imaginaire de la transformée de Fourier de cette fonction (figure IV.6) sont respectivement proportionnelles à la composante d'absorption et à la composante de dispersion du spectre RMN de l'objet en question. Ce résultat est, en effet, en parfait accord avec celui des travaux de Blümich et Ziessow [Blümich, 1982].
  97 ~(ml)
  0.08
  N=64
  0.07 0.06
    
  0.05
  
  0.04
         
  0.03 0.02 0.01
  o '"
  -0.01 '--_ _--'-_ _ _-'-_ _ _'--_ _--'-_ _ _1-_ _- 1_ _ _--' o 10 20 30 40 50 60 70
  :.
  10 8
  . "
  . '
  ,, ' :
  ,';
  "
  ",
   ;  6
    
  4
   
   .
  ~
  :
  ;
    .~ . :
  2
  :'
  '
  ;
  ,
  . . 
  
  ".:
  . ï  " .:
  :
  o ".:
  -2
  80 70
  o
  0
  Figure IV.5 Les coefficients du premier (filtre linéaire) et le deuxième (filtre quadratique) noyaux de Volterra pour une mémoire de système N=64.
  98 o
  -0.05
  Absorption -0.1
  -0.15
  -02
  -0.25
  o
  - 10
  5
  10
  Fréquence en kHz
  0.03
  Dispersion 0.02
  0.01
  0
  -0.01
  -0.02
  -0.03
  -O.ot -10
  -5
  o
  5
  10
  Fréquence en kHz
  Figure IV.6 La partie réelle (absorption) et la partie imaginaire (dispersion) de la transformée de Fourier du premier noyau de Volterra.
  99 111.2.3.2. Discussions Nous avons donc montré qu'un système RMN peut être modélisé par la mise en parallèle d'un filtre linéaire et d'un filtre quadratique. Actuellement, le traitement et l'analyse du signal RMN considèrent uniquement la contribution du filtre linéaire. Le filtre quadratique, qui est un complément de la partie linéaire, permettrait d'optimiser ce traitement et apporterait des informations complémentaires contenues dans un signal RMN. Pour avoir une idée de l'influence du filtre quadratique dans le modèle établi, nous avons évalué les contributions des deux noyaux de Volterra dans la réponse totale du système. La figure IV.7 montre que la contribution du second noyaux (filtre quadratique) est très faible comparée à celle du filtre linéaire. Ce résultat est prévisible du fait qu'en résonance magnétique, la réponse du système change de signe avec le changement du signe de l'excitation. La réponse est donc une fonction impaire de l'excitation et tous les termes pairs de la série de Volterra doivent être théoriquement nuls. En pratique, tronquer la série de Volterra à l'ordre 2 conduit à une contribution faible non nulle du filtre quadratique.
  n
  Figure IV.7 Les contributions du filtre linéaire (en pointillé) et du filtre quadratique (++) dans la réponse totale (traits continus) du système RMN.
  100 Pour améliorer ce modèle, il conviendra donc de faire appel au noyau de Volterra de 3ème ordre.
  111.3. Identification par série de Volterra de 3ème ordre Le développement à l'ordre 3 de la série de Volterra conduit à l'équation suivante:
  N-I
  d(n)
  =
  N-I
  Lhlml)x(n - ml) + L N-I
  +L
  N-I
  Lhz{m p m 2 ) x(n - ml)x(n - ml) +
  (IV.23)
  N-I
  N-I
  L
  Lhlm p m 2 , m 3 ) x(n-m l )x(n-m 2 ) x(n-m 3 )
  IIII~O m2~O m3~O
  Dans ce cas, le schéma général d'identification est celui présenté dans figure (IV.8). Il s'agit maintenant d'identifier les trois premiers noyaux de Volterra. Les équations de mise à jour des coefficients par l'algorithme de gradient sont alors [Asfour 2002] :
  e(n)
  = den) - H{ (n)X I (n) - H; (n)X 2 (n) - Hi (n)X 3 (n)
  °
  (IV.24)
  fIl (n + 1) = fIl (n) + 1(n)X I (n)e(n)
  (IV.25)
  fI 2 (n + 1) = fI 2(n) + 02 (n)X 2(n)e(n)
  (IV.26)
  fI3 (n + 1) = fI3 (n) + 03 (n)X 3(n)e(n)
  (IV.27)
  où X 3 (n) et H 3 (n) sont respectivement le vecteur d'entrée du filtre de 3ième ordre et le vecteur de ses coefficients. Ils/sont définis à l'instant n par :
  101 SystèmeRMN Signal de référence
  x(n)
  den)
  Filtre de Volterra signal d'entrée Filtre linéaire Filtre quadratique
  erreur e(n) = den) - d en) N-l
  d(n) = L: ~(~ :n)x(n - m1 )+ N- l
  + L:
  N- l
  L:
  N-l N-l
  L:
  L: h2 (~,m2 :n)x(n - m1)x (n -m2)
  N-l
  L: ~("'t,m2,m3: n)x(n-m1)x(n-m2)x(n-m3)
  ml ~ O ml ~O m3~ O
  Figure IV.8
  H3(n)
  Schéma général d'identification par filtre de Volterra de 3ième ordre.
  = [h3(0,0,0),h3(0,0,1),h3(0,0,2),  ,h3(N -1,N -1;N -1)]
  (IV.28)
  X 3 (n) = [x\n),x 2 (n)x(n -1),x 2 (n)x(n - 2),  ,x 3 (n - N + 1)]
  (IV.29)
  Par analogie avec les pas d'adaptation donnés par les équations IV.20 et IV.21, nous définissons le pas d'adaptation 8 3 (n) comme: 8 3 (n)
  =
  ~3
  3
  (NCY x (n»
  Ce modèle a été validé en reconsidérant l'objet ID décrit au chapitre 1 et ce pour une impulsion sine 90° de 1ms de durée comme signal d'entrée. La mémoire du système a été choisie à N=64 pour permettre un temps de calcul raisonnable et une convergence rapide de l'algorithme. Le résultat de l'identification est données dans la figure IV.9 où sont superposés le signal de référence et le signal de sortie du filtre de Volterra de 3ième ordre.
  102 0.8
  OB
  0.4
  02
  400
  450
  500
  550
  600
  650
  700
  n Figure IV.9 Le signal de référence (traits continus), et le signal à la sortie du filtre de Volterra pour une mémoire de système de 64.
  Le système RMN est donc équivalent à l'association en parallèle de trois filtres. La contribution de chacun de ces filtres dans la réponse totale est donnée dans la figure IV.l 0 où il est aisé de constater que le terme quadratique s'annule. La contribution maximale, en terme
  d'amplitude, du filtre de 3ième ordre est d'environ 6% dans la réponse totale confirmant les premiers résultats mis en évidence par Hoult [Hoult, 1979).
  103 0.9
  Partie linéaire
  OB 0 .7
  0.6
  2ème ordre
  0.5 o .~
  3ème ordre
  0.3 02 0.1
  ~oo
  ~20
  ~O
  ~60
  ~80
  500
  520
  5~
  n
  Figure IV.I0 Les contributions des trois filtres de Volterra dans la réponse totale du système RMN. En conclusion, la représentation par série de Volterra est très attractive dans les applications de filtrage étant donné que la sortie du système peut être considérée comme un produit de convolution entre les coefficients des noyaux et des combinaisons du signal d'entrée, soit un «filtrage linéaire ». L'inconvénient majeur de cette représentation réside dans le fait qu'un grand nombre de coefficients est généralement nécessaire pour décrire la non linéarité. Ceci a pour effet de rendre cette représentation peu commode pour les applications du traitement de signal. TI est alors préférable de s'orienter vers d'autres modèles qui seraient plus parcimonieux dans l'usage de leur coefficients. Nous introduisons donc ci-après un modèle récursif non linéaire qui pourrait constituer une alternative à la représentation par série de Volterra. Ce modèle sera détaillé dans le chapitre V et ne traitera que des signaux RMN réels acquis à l'aide de l'instrument que nous avons développé. Cependant, nous donnons pour permettre une comparaison avec le modèle de Volterra, la définition de ce modèle et les résultats obtenus sur des signaux RMN simulés.
  104 Modèle du filtre à réponse impulsionnelle infinie non linéaire
  IV.
  D'une façon générale, la relation entrée-sortie d'un filtre à réponse impulsionnelle infinie (RID non linéaire de mémoire N est décrite par une équation récursive de la fonne :
  M
  d(n)
  = L/~(d(n -l),d(n -
  2), d(n - N + l),x(n),x(n -l),  ,x(n - N + 1))
  (IV.30)
  i=1
  où p;(e,e,  ,e) est un polynôme d'ordre i des quantités entre parenthèses. Ce modèle du filtre RH non linéaire pennet de décrire les caractéristiques non linéaires d'un système avec moins de coefficients qu'un modèle basé sur la série de Volterra [Mathews, 1991] Le plus simple des filtres RH non linéaires est celui dont l'entrée et la sortie sont liées par l'équation :
  N-I
  d(n)
  =
  N-I
  N-I
  N-I
  LC(i)d(n-i) + L
  Lb(i,j)d(n- j)x(n-i)+ La(i)x(n-i)
  i=1
  j=1
  i=O
  (IV.3l)
  i=O
  où cri), b(i,j), a(i) sont les coefficients du modèle. Par analogie avec les appellations employées pour les coefficients d'un filtre (RII) linéaire, les coefficients cri) et a(i) seront respectivement appelés coefficients auto-régressifs et coefficients de moyenne mobile. Les coefficients b(i,j) sont les coefficients non linéaires. L'équation IV.31 traduit le fait que la sortie du système d(n) à un instant donné n est calculée à partir des N dernières valeurs de l'entrée et des N-l dernières valeurs de la sortie. L'identification du système consiste alors à estimer les coefficients cri), b(iJ) et a(i) du filtre. Ces coefficients sont donnés, pour l'objet ID décrit au chapitre I, dans la figure IV.ll. Le signal d'entrée étant une impulsion sinc 90° de lms de durée. L'algorithme de gradient, qui sera détaillé dans le chapitre suivant, a été utilisé pour l'estimation des coefficients. La mémoire du système a été fixé à N =24.
  105 a(i) 0.06
  +-
  0.04
  +.
  '* * * * *
  0.02 0
  31<-
  *
  '" *
  31<-
  +.
  * '* * * * *
  31<;jo
  31<-
  -0.02 -0.04
  * "* 0
  5
  10
  15
  ;
  20
  25
  c(i) 0.2
  *
  0.15 0.1
  '*
  '*
  '*
  +
  *'
  0.05 0 -0. 05 0
  5
  t.-
  +.
  *
  31<-
  .;<-
  '"
  10
  +-
  '*
  t.-
  +.
  * *' *
  "* "" 20
  15
  +-
  *" 25
  i
  6
  5 4 3
  2
  o -1 -2
  25 25
  o
  0
  ;
  Figure IV.II Les coefficients du filtre RIT non linéaire obtenus l'algorithme de gradient. La mémoire du système est N=24.
  106 Le signal de sortie du filtre RIT ainsi calculé est donné en même temps que le signal de référence dans la figure N .12.
  OB
  0.6
  0 .4
  02
  ~2--------~------~------~~------~------~-------400 450 500 700 550 600 650
  Figure IV.12 Le signal de référence (traits continus) et le signal de sortie du filtre RIT non linéaire (en pointillés) pour N=24. On peut dire maintenant que le système RMN est équivalent, au sens de l'erreur quadratique moyenne, à un filtre RIT non linéaire [Asfour-JMR, 2000], [Asfour, 2000].
  Al' issu de chapitre, nous pouvons conclure que la structure RIT non linéaire permet de modéliser la non linéarité de systèmes RMN avec un nombre moins élevé de coefficients que la représentation par série de Volterra. C'est précisément ce modèle qui sera désormais retenu dans l'identification de systèmes RMN réels. TI sera donc exposé en détails dans le chapitre suivant.
  107 Chapitre V
  Modélisation des systèmes RMN par filtre à réponse impulsionnelle infinie (RII) non linéaire: Validation sur des signaux réels
  108 Introduction Nous avons montré dans le chapitre précédent, que la représentation de systèmes RMN par le modèle de série de Volterra n'est pas appropriée pour décrire la non-linéarité en raison du grand nombre de coefficients à estimer. Nous avons donc introduit le modèle du filtre RII non-linéaire et démontré qu'il permet les mêmes performances et ce avec un nombre moins élevé de coefficients. Nous allons donc, dans ce chapitre, détailler le modèle du filtre RII non linéaire. Nous verrons tout d'abord les modifications apportées pour réduire le nombre de coefficients. Nous présentons ensuite l'algorithme d'adaptation mis en oeuvre et les approches du filtrage adaptatif RII associées à savoir l' approche d'erreur d'équation et l'approche d' erreur de sortie . Le modèle sera finalement validé sur des signaux RMN réels acquis à l'aide de l'instrument développé et décrit précédemment (voir chapitres II et III).
  Il
  Modèle du filtre RII non linéaire et méthode d'identification Il.1
  Description du modèle Le modèle du filtre RII proposé est une variante du modèle décrit par l'équation 11.31.
  Elle est réalisée en agissant sur la mémoire du système, préalablement notée N, d'une manière à ce que le nombre des coefficients de chaque partie du filtre puisse être choisi indépendamment les uns des autres. L'équation décrivant ce modèle de filtre est ainsi de la forme [Asfour, 2001] :
  K-I
  L'-I K'-I
  den) = LC(i) den - i) + L ;=1
  ;=0
  L-I
  Lb(i,j) den - j) x(n - i) + La(i) x(n - i)
  (V.1)
  ;=0
  j=1
  où L, L' , K, et K' sont des nombre entiers positifs. L'intérêt de ce modèle réside dans le fait qu'il permet d'approcher une large classe de systèmes non linéaires et ce avec un nombre fini de coefficients. De plus, la plupart des concepts développés le concernant peuvent être facilement étendus à des modèles non linéaires plus généraux.
  Il.2
  Les approches du filtrage adaptatif RII non linéaire
  109 L'identification des systèmes RMN par filtrage adaptatifRII non linéaire est basée sur le schéma général de la figure V.l. Le signal d'entrée x(n) est appliqué au système RMN considéré qui produit en sortie le signal de référence d'(n) .
  Bruit de mesure x(n)
  SystèmeRMN
  d(n)
  d'(n)
  + e(n)
  ~
  +'"~
  retard
  1
  "
  d(n)
  Filtre RII non linéaire "
  y(n) = d(n)
  1
  ou
  d(n) K -I
  d(n)
  = L: ê(i: n) y (n - i) + '~ I
  Figure V.I linéaire.
  L' - I K'-I
  L- I
  1= 0
  '~ O
  L: 2: b (i, j:n) y(n - j)x(n - i) + L: â(i: n) x(n-i) }=I
  Schéma général pour l'identification des système RMN par filtrage adaptatif RIT non
  Pour le traitement des signaux réels, il convient, bien entendu, de prendre en compte le bruit de mesure qui s'additionne au signal de référence d'(n) pour donner le signal d(n) qui est réellement utilisé dans le schéma d'identification. Dans la suite, nous considérons que le niveau de ce bruit est suffisamment faible pour ne pas influencer les performances de l'algorithme adaptatif. Le filtre adaptatif va donc essayer d'estimer le signal de référence d(n). Cette estimation est donnée à chaque instant n par :
  K-I
  d(n) = Lê(i : n) y(n - i) + ;=1
  L'-I K'-I
  L-I
  ;=0
  ;=0
  L Lb(i,j: n) y(n - j) x(n - i) + Lâ(i: n) x(n - i) j=1
  110
  (V.2) où ê(i: n), â(i : n) et b(i, j : n) sont les coefficients du filtre adaptatif à l'instant n. Ces coefficients sont itérativement mis à jour pour minimiser l'erreur quadratique moyenne. Le calcul itératif de ces coefficients peut être fait selon deux approches distinctes. Dans la première approche dite d'erreur d'équation, les échantillons x(n) et d(n) du signal d'entrée et du signal de sortie respectivement sont utilisés pour estimer la sortie den) du filtre adaptatif. Ceci correspond, évidemment, à remplacer dans l'équation V.2 le signal y(n) par
  d(n). L'approche d'erreur de sortie considère, quant à elle, les échantillons de x(n) et les échantillons passés de den) pour calculer l'estimation actuelle d(n) , soit donc en utilisant l'équation V.2 avec y(n)= d(n) . Les principaux mérites de ces deux approches seront discutés plus loin dans ce chapitre. Pour plus de détail les concernant, le lecteur pourra se rapporter à la littérature [Shynk, 1989].
  Il.3
  L'algorithme d'adaptation Il.3.1 Les équations de mise à jour
  Nous avons employé, comme pour l'identification par série de Volterra, l'algorithme de gradient pour identifier les coefficients du modèle. Les équations de mise à jour des coefficients ê(i : n), â(i: n) et b(i, j : n) du filtre sont:
  â(i: n+l) =â(i: n) +ba x(n+l - i)e(n+l)
  (V. 3)
  b(i, j: n + 1) = b(i,j: n) + bby(n + 1- j)x(n + 1- i)e(n + 1)
  (V.4)
  ê(i: n + 1) = ê(i, n) + bcy(n + 1- i)e(n + 1)
  (V.5)
  où 0a' Oc et Ob désignent les pas d'adaptation. En faisant appel aux notations matricielles l , l'expression de den) donnée dans l'équation V.2 peut s'écrire par :
  "AT
  T'" T
  ,
  den) = A (n)X L (n) + trace[B(n)(YK • (n)X L• (n))] + C (n)YK (n)
  111
  (V.6) Dans cette expression: Â(n) et ê(n) sont respectivement les vecteurs des coefficients moyenne mobile et auto-
  régressif à l'instant n :
  Â(n)
  = [â(O: n),â(l: n),  ,â(L -1: n)f
  ê(n)
  = [ê(l: n),ê(2: n), ,ê(K -
  (V.7) (V.8)
  l;nf
  B(n) représente la matrice, de dimension L' x (K' - 1) , des coefficients non linéaires à
  l'instant n :
  ~
  b(O,l : n) ~
  b(l,l: n) Ê(n)
  ~
  b(O,K' - 1: n)
  b(O,l : n) ~
  b(l, K' - 1: n)
  b(1,2 : n)
  =
  (V.9)
  b(L ' - 1,1: n )
  b(L ' -1,2: n )
  b(L ' -l,K' -l:n)
  X L(n ),XL,(n),YK (n) et YK, (n) sont les vecteurs d'entrée du filtre adaptatif à l'instant n. Il
  sont définis par :
  XL (n) = [x(n), x(n -1),  , x(n - L + l)y
  (V. 10)
  = [x(n),x(n -l),  ,x(n -
  (V.l1)
  XL' (n)
  L' + l)y
  YK(n)
  = [y(n -l),y(n -
  2),  ,y(n - K + l)y
  (V.12)
  YK,(n)
  = [y(n -l),y(n -
  2),  ,y(n - K' + l)y
  (V. 13)
  Avec ces nouvelles notations, les équations de mise àjour s'écrivent alors:
  e(n + 1) =d(n + 1) - ÂT (n)Xdn + 1) - trace[Ê(n)(YK,(n + l)XI,(n + l))]_ê T (n)YK(n + 1) Â(n + 1) = Â(n) + baXJn + 1)e(n + 1)
  (V.l4) (V.l5)
  1 (f et trace désigne respectivement la transposée d'une matrice et la somme des éléments diagonaux d'une matrice carrée.
  112 ê(n + 1)
  = ê(n) + <5 c Y K (n + l)e(n + 1)
  Ê(n + 1) = Ê(n) + <5 b X L ,(n + l)YK .(n + l)e(n + 1)
  (V.16) (V.17)
  Lorsque les nombres L, L', K' et K sont fixés, la stabilité de l'algorithme et la qualité des résultats ne dépendent que du choix des pas d'adaptation.
  Il. 3.2 Conditions de stabilité et pas d'adaptation variables Nous rappelons préalablement que les conditions de stabilité dans une réalisation du filtrage adaptatifRII linéaire utilisant l'algorithme de gradient sont données par les inégalités suivantes [Bellanger, 1989] :
  2
  0< Sa < - -- -- - - (L - 1)0-~ + (K - 1)0-;
  (V.1 8)
  2 0 < Sc < - -- - - - - (L - 1)0-~ + (K - 1)0-;
  (V.19)
  où o-~ et 0-; sont respectivement les puissances des signaux x(n) et y(n) En revanche, les conditions de stabilité de l'algorithme d'adaptation explicité
  Cl-
  dessus doivent être établies car nous traitons de filtrage adaptatif non linéaire.
  Il.3.2.1
  Conditions de stabilité sur les pas d'adaptation
  L'examen de l'équation V.14 montre que la suite d'erreurs e(n) est calculée en utilisant les coefficients du filtre avant leur mis à jour. Elle est, par conséquent, appelée suite des erreurs « a priori ». L'erreur dite « a posteriori» définie par l'équation:
  &(n + 1) = d(n + 1) -
  ÂT
  (n + l)X L (n + 1) - trace[Ê(n + 1)(YK , (n + l)X[.(n + 1))] - ê T (n + l)YK (n + 1) (V.20)
  est, en revanche, calculée à partir des coefficients mis à jour, i.e. en utilisant les valeurs des coefficients obtenues une fois terminées les opérations de mise à jour décrites par les équations (V.15), (V.16) et (V.17). La relation entre ces deux types d'erreur peut s'écrire comme:
  113 e(n + 1) - s(n + 1)
  (V.21) En prenant en compte les équations V.1S, V.16 et V.17, il est possible d'exprimer l'erreur a posteriori en fonction de Gelle a priori par :
  e(n+1)=e(n+1)
  {1 -t5a X[(n+l)XJn+1)
  - t5 b trace[(X L.(n+1)yJ,(n+1))(YK,(n+1)XI.(n+1))]
  -t5c Y: (n + l)YK (n + 1)} (V.22)
  L'algoritlune d'adaptation peut être considéré stable si l'espérance mathématique de la valeur absolue de l'erreur a posteriori est plus petite que celle de l'erreur a priori, ce qui est logique puisque s(n + 1) exploite davantage d'information [Bellanger, 1989]. Pour exprimer cette condition de stabilité sous une forme mathématique (ou inégalité), nous supposons que la partie linéaire du filtre agit indépendamment de sa partie non linéaire. Il est donc possible de considérer que la contribution dans l'erreur totale, qu'elle soit a priori ou a posteriori, de la partie linéaire est indépendante de la contribution de la partie non linéaire. Si nous désignons maintenant par sL(n) et SNL(n)
  (respectivement eL(n) et
  e NL (n)) les contributions de la partie linéaire et de la partie non linéaire dans l'erreur a posteriori (respectivement a priori), nous obtenons les deux conditions de stabilité que les pas
  d'adaptation doivent satisfaire, soit:
  (V.23) et (V.24)
  La première condition (V.23) permet de définir l'intervalle des valeurs permises pour les pas 8 a et 8c ' et d'établir les conditions de stabilité données par les inégalités V.l8 et V.19. La deuxième condition (V.24) permet de calculer les limites à respecter par le pas d'adaptation non linéaire t5b • Pour cela, nous la réécrivons en considérant l' équation (V.22) :
  114 E[leNJn+l)-eNJn+l) 1] = jl-(Vrace{E[(Xdn+1)Y;'(n+l))(YK,(n+l)X[.(n+l))]}j < 1 (V.25) Nous faisons ensuite l'approximation suivante:
  (V.26)
  où
  r xy (0)
  Pour
  des
  est la fonction d'intercorrelation des signaux x(n) et y(n) à l'instant d' origine n=O. signaux
  à
  r XY (0) ~ min(L',K')a;y
  bande
  limitée,
  cette
  fonction
  peut
  être
  écrite
  comme :
  où a ;y désigne la puissance moyenne. Il en résulte la condition de
  stabilité suivante:
  (V.27) soit encore:
  2 0 < 8b < - - - - - - min(L ', K ')( a ;y)2
  (V.28)
  En pratique, il serait prudent de prendre une marge de sécurité correspondant à un facteur de quelques unités pour pouvoir prendre en compte des approximations faites dans le calcul.
  Il.3.2.2
  Les pas d'adaptation variables pour les signaux RMN
  Pour permettre le traitement de signaux non stationnaires, comme c'est le cas des signaux RMN, les algorithmes adaptatifs doivent être capables de suivre l'évolution des paramètres de ces signaux. Dans le contexte actuel de l'identification de système RMN, les deux signaux d'entrée du filtre RH sont non stationnaires. L'apparition d'instabilités est, par conséquent, inévitable si des valeurs constantes sont utilisées pour les pas d'adaptation. Pour remédier à ce problème, nous avons introduit, comme dans le cas de l'identification par série de Volterra, des pas d'adaptation variables en fonction de n :
  115 8 a(n)
  où
  =
  Pa
  (V.29)
  (L -1)a; (n) + (K -1)a; (n)
  8c(n) = Pc (L-l)a;(n)+(K -1)a;(n)
  (V.3D)
  8b(n) = Pb min(L',K')(a;y (n))2
  (V.31)
  Pa' Pc
  et
  Pb
  sont des constantes positives ; (}; (n) et (}~(n) sont respectivement des
  estimations des puissances des signaux x(n) et y (n) , et a;y (n) est une estimation de la puissance moyenne à l' instant n. Ces variables sont définies par les relations suivantes:
  (V.32)
  (V.33)
  (V.34)
  où a;o a;o et a;yO sont des constantes positives permettant d'éviter une division par zéro. Ces estimations de puissance sont calculées sur des fenêtres temporelles de largeur respective
  III
  Validation du modèle sur des signaux réels RMN 111.1
  Exemple de signaux réels
  La figure V.2 montre un exemple de signaux RMN acquis à l'aide de l'instrumentation développé.
  116 -0:2
  -0.4
  0 ,-* --,--4-'-°--,.--,-,6-'-°-----,-l.80-. --~10'-:.0---12-'-0----:-*.,.,,1L40,.-*, ,-J ,0.8 L-,.--2'-:.
  (n) Figure V.2 Exemple de signaux RMN réels. L'excitation RF est une impulsion rectangulaire 90° de 2.5ms de durée et le signal FID est observé pendant une période d'environ Ilms.
  Ce signal, de type Off-resonance, est obtenu pour un tube cylindrique rempli d'eau déminéralisée. Le champ magnétique statique est réglé à 0.11 T, ce qui correspond à une fréquence de Larmor d'environ 4.783MHz. L'excitation RF est une impulsion 90° de forme rectangulaire et de 2.5ms de durée. La fréquence d'échantillonnage dans la bande de base est de 10kHz et la durée de l'observation du signal est de 11ms, soit donc un nombre total de 110 échantillons acquis sur cette durée.
  111.2
  Résultats d'identification
  Les paramètres L, K, L' et K' doivent être choisis pour permettre une meilleure description du système RMN. Pour cela, nous nous somme fixés les critères suivants. Il est d'abord important de garder à l'esprit que le nombre de coefficients doit être raisonnable pour réduire le temps du calcul et simplifier l'exploitation du modèle dans des applications éventuelles. Contrairement au modèle de Volterra, nous choisissons un nombre de coefficients non linéaires moins élevé que celui des coefficients linéaires pour traduire le
  117 fait que contribution de la partie non linéaire dans la réponse totale du système est moins importante comparée à celle de la partie linéaire. Ces coefficients doivent, bien entendu, être choisis de manière à minimiser l'erreur quadratique moyenne entre le signal de référence et la sortie du filtre. Toutes ces considération nous ont conduit au choix des valeurs suivantes : L=12, K= 22, L'
  =8
  et K'
  = 2,
  soit donc 41 coefficients à estimer. Dans l'exemple du signal RMN
  donné dans la figure V.2 où nous disposons d'environ 110 échantillons du signal, l'estimation des coefficients du filtre correspondant est significative. Cependant, la qualité du modèle qui est fonction de la redondance nombre de mesures/nombre de paramètres peut être améliorée en terme de précision d' estimation des coefficients en faisant appel à des mesures complémentaires du même signal RMN. Ainsi, dans la mise en oeuvre pratique de l'algorithme d' adaptation, nous n'utilisons pas qu'une seule excitation RF et le signal RMN résultant mais les signaux correspondant à la répétition de cette excitation. L' impulsion d' excitation est ainsi répétée périodiquement et le signal RMN est acquis à la fin de chaque impulsion. La fréquence de répétition devra être choisie adéquatement pour permettre au système d' atteindre son état d'équilibre avant l' application de l'impulsion suivante. A chaque impulsion, l'algorithme adaptatif effectue la mise à jour des coefficients de filtre en prenant pour valeurs d'initialisation les valeurs de ces même coefficients à fin de l' impulsion précédente. A la première impulsion, ces coefficients utilisent des valeurs de départs fixées par l'utilisateur. Nous montrerons par la suite que les valeurs des coefficients se stabiliseront après un nombre donné de répétitions de l'impulsion d'excitation [Asfour, 2001]. Examinons dans ces conditions l'estimation de ces coefficients selon l'approche d'adaptation basée sur l'erreur d'équation. Ceux-ci étant initialisés à zéro et les valeurs de différents paramètres d'adaptation, conduisant au minimum de l'erreur quadratique moyenne entre le signal de référence et le signal à la sortie du filtre, sont: fle
  = 0.265,
  Nx
  = 5,
  Ny
  = 10,
  N xy
  = 5,
  (J'~o
  = 0.030,
  (J'~o
  = 0.058,
  fla
  = 0.680,
  (J'~yO
  = 0.085
  flb
  = 0.010,
  et ~ = 3.
  Le signal de sortie du filtre résultant est comparé au signal de référence dans la figure V.3 qui confirme, d'ailleurs, que le système RMN est équivalent à un filtre RH non linéaire au sens de l'erreur quadratique moyenne. Notons que le signal de référence considéré dans cette analyse n'est autre que le signal RMN réel donné dans la figure V.2 et que les valeurs initiales de la sortie du filtre, i.e. pendant la durée de l'impulsion d'excitation, ont été fixées à zéro.
  118 0.8 0.6 0.4 0.2
  -0.2 -0.4 -0.6 -0. 8 ' -_ _- L_ _ _-L-_ _ _.L.-_ _-----''-_ _--'-_ _ _---'-_ _ _- ' - - - - - ' 140 20 40 60 80 100 120
  n Figure V.3 Le signal de référence (traits continus) et le signal à la sortie du filtre Rn non linéaire (en pointillés). Les coefficients du filtre sont estimés à l'aide de l'algorithme d'adaptation selon l'approche d'erreur d'équation.
  Nous devons préciser que ce résultat a été obtenu pour 30 répétitions du signal d'excitation qui consiste, nous le rappelons, d'une impulsion 90° de forme rectangulaire et de 2.5ms de durée. Pour avoir une idée de l'influence du nombre de répétions de l'impulsion RF, soit donc de l'influence du nombre de mesures du même signal RMN, sur la qualité du modèle, nous présentons dans la figure V.4 les variations de l'erreur quadratique moyenne (EQM) en fonction du nombre de répétitions de l'impulsion d'excitation. Cette erreur peut être considérée constante au delà de 30 répétitions de l'impulsion RF et indique ainsi qu'au moins 30 mesures du même signal RMN seront nécessaires pour obtenir l'estimation optimale, en terme de précision, des coefficients du filtre. Malgré ce nombre assez élevé, le temps de calcul reste raisonnable et autorise l'implantation de l'algorithme d'adaptation en temps réel.
  119 X 10.3
  EQM 4.5
  4
  3.5
  3
  2.5
  2
  1.5
   . . . . . . . . . . . . . . . . . . . . . . . 5
  10
  20
  15
  25
  30
  35
  40
  45
  50
  Numéro de /'impulsion
  Figure V.4 L'erreur quadratique moyenne (EQM) en fonction du nombre de répétitions de l'impulsion RF.
  L'approche d'erreur d'équation présente l'avantage de l'unicité de la solution du problème, c'est-à-dire que l'algorithme d'adaptation converge vers le seul minimum global de la fonction d'erreur à minimiser. Son intérêt est, cependant, limité aux cas où le bruit additif au signal de référence est de faible niveau. Dans le cas contraire, cette approche conduit généralement à une estimation biaisée des coefficient du filtre [Mathews, 1991], [Shynk, 1989] et [Michaut, 1992]. Théoriquement, il est alors préférable d'employer l'approche d'erreur de sortie qui est l'approche correcte. Cette approche qui utilise uniquement, comme nous l'avons vu, les échantillons passés de
  J (n ) pour
  estimer l'échantillon présent (actuel) de
  J (n)
  devra
  permettre l'obtention d'une estimation non biaisée des coefficients. Pour tester cette approche, les valeurs initiales de ces coefficients s'avèrent être les paramètres les plus cruciaux de l'algorithme. En fait, s'ils ne sont pas proprement initialisés, l'algorithme d'adaptation peut bel et bien converger vers un des minima locaux que possède
  120 la fonction d'erreur à minimiser. Ce minima local pouvant être, bien entendu, très différent du minimum global. Les valeurs initiales des coefficients qui garantissent la convergence de l'algorithme vers le minimum global peuvent être fixées soit empiriquement soit en utilisant des valeurs proches de celles des coefficients issus de l'algorithme d'adaptation basé sur l'approche d' erreur d'équation. La figure V.S montre le résultat de l'identification basée sur l'approche d'erreur de sortie en se basant sur cette deuxième méthode pour l'initialisation des coefficients
  0.8 0.6 0.4 0.2
  -0.2 -0.4 -0.6 -0.8 L -_ _----L-_ _ _L-_ _---1_ _ _---L-_ _ _- ' - -_ _---'-_ _ _- - - L - - ' 100 120 140 20 40 60 80
  n Figure V.S Le signal de référence(traits continus) et le signal à la sortie du filtre RII non linéaire (en pointillés). Les coefficients du filtre sont estimés par l'algorithme d'adaptation basé sur l'approche d'erreur de sortie.
  En conclusion à ce chapitre, nous avons développé un modèle pour l'identification de systèmes RMN et nous l'avons validé sur des signaux réels. Ce modèle nous a permis d'approcher un système RMN par un filtre Rn non linéaire. Malgré les difficultés non entièrement résolues à ce jour dans la mise en oeuvre de ces filtres, notamment les problèmes de stabilité, nous avons pu conduire la mise en place d'un
  121 algorithme d'adaptation. Cet algorithme pourra constituer une contribution . dans la problématique générale de l'identification des systèmes non linéaires et dans l'utilisation, encore peu développée, de structures RH non linéaires dans le domaine du traitement de signal. Nous signalons aussi que notre approche d'identification, a été conduit uniquement dans le cas d'une expérience de base en RMN où le terme non linéaire, comme nous avons pu le constaté, ne contribue que d'environ de 6% dans la réponse totale du système. Notre modèle doit être comparé avec les modèles paramétriques linéaires, actuellement en utilisation, pour pouvoir apprécier ses avantages. Ajoutons finalement que ce traitement non linéaire pourra trouver tout son intérêt dans des expériences plus complexes de la RMN telles que celles avec des impulsions multiples ou encore celles de la résonance stochastique ou les effets de la non linéarité sont plus importants.
  122 Conclusions et perspectives
  123 Les travaux de recherche et de développement, menés durant cette thèse, sur la mise en place d'un imageur numérique, ont été fructueux et ont donné naissance à un produit quasi fini . Les performances de l'imageur développé, en particulier celles du récepteur et du synthétiseur numériques, constituent une avancée importante dans le développement instrumental et le traitement du signal pour l'IRM dédiée à bas champ. Outre son coût peu élevé et le bon RSB atteint, le caractère peu encombrant de l'électronique d'acquisition et de traitement que nous avons élaboré ouvre la voie à une intégration facile dans un cabinet médical ou bien sur les lignes de production pour le contrôle de qualité agroalimentaire. Les images esquissées obtenues par la séquence d'écho de gradient peuvent être améliorées. L'imageur possède le potentiel d'intégrer une large variété de séquences d'imagerie. Pour pouvoir apprécier la qualité des images en diagnostic médical ou en contrôle agroalimentaire, il y a bien évidemment lieu de les comparer à celles d' un imageur bas champ utilisant une électronique de réception analogique. La solution numérique de réception du signal RMN laisse envisager des perspectives de recherche encore plus poussées pour rehausser davantage le RSB. L'application, par exemple, des techniques de filtrage adaptatif multivoies au traitement du signal RMN permettrait d'apporter des améliorations supplémentaires sur le RSB pour éviter le moyennage du signal et réduire ainsi le temps d'imagerie. Concernant la deuxième partie du travail présenté dans ce manuscrit, nous avons proposé et comparé deux modèles pour l'identification des systèmes RMN. Ces modèles trouvent leur intérêt dans le fait qu'ils permettent de prendre en considération les caractéristiques non linéaires, qui sont jusqu'à présent négligées, de tels systèmes. Ainsi, nous avons montré qu'il était possible de modéliser un système RMN aussi bien par un filtre de Volterra de troisième ordre que par un filtre non linéaire à réponse impulsionnelle infinie (RIn. Néanmoins, nous avons retenu ce dernier grâce à son potentiel de décrire la non linéarité des systèmes RMN avec moins de coefficients que le modèle de Volterra. L'apport du modèle du filtre RII dans l'amélioration du traitement et de l'analyse du signal RMN reste à étudier. Pour cela, il est nécessaire d'entreprendre une étude de comparaisons avec les méthodes paramétriques linéaires de traitement du signal RMN qui sont actuellement utilisées et qui pourront d'ailleurs trouver leur limitation dans les expériences de la RMN stochastique ou encore dans celles à impulsions multiples. D'autres études peuvent être également envisagées pour explorer la possibilité de relier les coefficients du modèle non linéaire aux propriétés physiques de l' échantillon
  124 modélisé. Ceci nécessite, bien évidemment, l'analyse approfondie du comportement des différents coefficients du modèle pour un nombre significatif d'échantillons.
  125 Annexe
  Le processeur ADSP21 060 SHRAC : Architecture et fonctionnement
  126 A.1
  Description générale Les processeurs ADSP-21060 SHARC d'Analog Deviees sont des DSP 32 bits virgule
  flottante hautes performances optimisés pour le traitement numérique du signal en temps réel. Ils sont aussi largement utilisés pour d'autres applications nécessitant une grande vitesse de calcul. La réputation de ce DSP doit principalement à sa capacité d' effectuer les calculs en virgule flottante. Il répond ainsi à des critères de qualités plus élevés que son équivalent en virgule fixe dans les domaines de : •
  Précision
  •
  Plage dynamique
  •
  Facilité et flexibilité d' utilisation
  Un schéma simplifié de l'architecture du processeur ADSP-21 060 SHARC est donné par la figure A.I. On peut constater sur ce schéma le niveau élevé d'intégration de ce DSP qui combine en un seul circuit intégré le coeur du processeur avec d'autres utilités. Celles-ci incluent en premier lieu la mémoire SRAM (Static Random Access Memory) qui peut être configurée par l'utilisateur en deux types de mémoire: une appelé « mémoire programme» et l'autre appelé «mémoire de donnée ». Cette mémoire peut être accédée via le port «processeur» par le coeur du processeur grâce aux quatre bus internes de celui-ci. •
  Les deux bus de programme PM : - PMD (program Memory Data) 48 bits : il est utilisé pour le transfert des instructions ou des données situées dans la «mémoire programme» entre celle-ci et le coeur du processeur et inversement. - PMA (Pro gram Memory Address) 24 bit: il est chargé du transfert des adresses des instructions et des données de la mémoire programme.
  •
  Le deux bus de données DM: - DMD (Data Memory Data) 32 bits: Il assure le transfert des données entre la mémoire de données et le coeur du processeur - DMA (Data Memory Address) 32 bits : il est utilisé pour transférer les adresses des données dans la mémoire de données
  127 En deuxième lieu, on remarque la présence d'un processeur indépendant IIO gérant les opérations de transfert de données entre la mémoire et les ports (ports séries, port externe, le contrôleur DMA) de communication avec des dispositifs extérieurs. Ce processeur peut accéder à la mémoire via le port « IIO » et ceci grâce aux bus IIO (IOA pour les adresses et IOD pour les données).
  CORE PROCESSOR
  , - - - DUAL-PORTED SRAM ------ TWQ INDEPENDENT DUAL*PORTED BLOCKS
  JTAG TEST & EMULATION
  32 DM ADDRESS BUS
  48
  CONTROL, STATUS & DATA BUFFERS
  Figure A.1
  Architecture du processeur ADSP21060 d'Analog Deviees.
  A.2. La mémoire SRAM Il s'agit donc de mémoire IMbits (1 mégabits) interne au processeur. Son architecture autorise le stockage d'une donnée ou d'une instruction dans la même location «physique », une particularité qui n'est pas forcément commune à tous les DSP, et qui a son grand intérêt pour la rapidité d'exécution de programme. En effet, cette mémoire est organisée en deux blocs. Chaque bloc peut stocker des combinaisons de code de programme (instructions sous format de 48 bits) et de données sur 32 bits. Cependant, l'accès à la mémoire est rendu plus efficace si l'un de ces blocs, que nous avons appelé «mémoire de données », est dédié au stockage de données -et qui utilise le bus DMD pour le transfert- alors que l'autre bloc, «appelé mémoire de programme» est chargée de stocke des instructions et des données en utilisant le bus PMD pour le transfert.
  128 Cette manière d'utilisation des bus de transfert (c'est-à-dire que chaque bus est dédié à un seul bloc mémoire garantit l'exécution d'une opération aux même temps qu'un transfert de deux données en un seul cycle d'horloge, à condition que l'instruction à exécuter soit disponible dans la mémoire « cache» (Instruction cache sur la figure A.l).
  A.3. Le coeur du processeur Le coeur du processeur contient trois unités de calculs, deux générateurs d'adresses et un séquenceur de programme avec mémoire cache.
  a) Les unités de calculs
  Ils sont en nombre de trois unités : - Une unité arithmétique et logique (UAL) : elle accomplit essentiellement les opérations arithmétiques et logiques standards en virgule fixe et en virgule flottante. - Un multiplieur accumulateur: il est chargé des opérations de multiplication en virgule fixe et en virgule flottante, ainsi que les opérations de multiplication/addition et multiplication/soustraction en virgule fixe. - Une unité de décalage: des opérations de manipulation de bits (décalage logique et arithmétique  ) peuvent être effectuées grâce à cette unité.
  Ces unités de calculs sont indépendantes et fonctionnent en parallèle. Chaque unité exécute une instruction de programme en un cycle d'horloge, et la sortie d'une unité peut être une entrée pour une autre dans le cycle suivant. Un ensemble de registres de données (Data Register File) de calculs en nombre de 16 permet d'assurer l'interface entre ces unités de calculs et les bus internes de données. Ces registres permettent aussi de stocker les résultats de calculs intermédiaires.
  b) Le séquenceur de programme et les générateurs d'adresses (DAGl, DAG2)
  Les générateurs d'adresses sont dédiés à fournir sur deux bus d'adresses PMA et DMA, les adresses des données et des instructions de programme situées dans la mémoire. En association avec ces générateurs, le puissant séquenceur de programme permet l' accès aux mémoires internes du processeur ou aux mémoires externes de la carte SRT-SHARC. Les
  129 deux générateurs d'adresses permettent au processeur de fournir simultanément les deux données nécessaires aux opérations à deux opérandes, et ceci pour un maximum d'efficacité et de puissance de calculs de 160MOPS (160 millions d'opérations par seconde). Du fait de son architecture, le processeur ADSP21 061 SHARC possède un haut degré de parallélisme spécifiquement adapté au traitement numérique de signaux et d'images. Au cours d'un seul cycle d'horloge, ce processeur peut effectuer les opération suivantes :
  •
  Générer les prochaines adresses du programme.
  •
  Aller chercher l'instruction suivante.
  •
  Effectuer un ou deux transferts de données
  •
  Exécuter une opération.
  •
  Recevoir ou transmettre des données à travers les ports séries
  L' exécution du programme s' effectue de manière linéaire, c'est-à-dire que les instructions sont donc exécutées par le processeur dans l' ordre selon lequel elles sont placées dans la mémoire de programme. Cependant, quelques fonctions permettent de rompre ce flous de programme. Celles-ci sont par exemple des fonctions de type boucles, sous-routines ou interruptions. Les interruptions, en particulier, sont causées par une variété de conditions internes ou externes au processeur. Les interruptions extérieures sont en nombre de trois. Elles sont causées par des dispositifs externes qui envoient une demande d'interruption sur l'une des entrées d'interruptions: IRQO, IRQ1 ou IRQ2. Parmi les interruptions internes au processeur, on peut citer les interruptions liées aux opérations de transfert du contrôleur DMA (Direct Memory Acess), et celles qui concernent la réception et la transmission des données par les ports série.
  A.4. Le processeur 1/0 et le port externe: Ce processeur fonctionne indépendamment du coeur du processeur. IL contient deux ports séries et un contrôleur DMA (Direct Memory Acces). Ce dernier constitue un moyen efficace et rapide pour effectuer le transfert d'un bloc de données entre la mémoire interne et une source de données extérieure (mémoire externe, PC, .  ). Les ports séries offrent un protocole de communication simple avec une large variété de périphériques de traitement numérique de signal. Les fonctions de réception de données
  130 par le DSP et de transmission vers la périphérique sont complètement indépendante pour un maximum de rapidité. Le transfert de données via le port série est dit de type synchrone. Autrement dit, le transfert des bits successifs d'un mot de données est rythmé à la cadence d'une horloge, il y a un transfert d'un bit à chaque cycle d'horloge. Cette horloge peut être générée par le DSP même, ou bien elle peut être fournie de l'extérieur. La fréquence de cette horloge peut aller jusqu' à la fréquence du processeur (40MHz). Ces ports séries sont employés dans notre système d'IRM pour assurer le transfert des données numériques (échantillon du signal RMN) du récepteur numérique vers la mémoire du DSP. Avec des bus parallèles, le port externe offre la possibilité d'interfacer le DSP avec d'autres processeurs (un autres DSP, PC,  ), avec une mémoire externe ou d' autres périphériques.
  A.5. Outils de programmation de la carte SRT-SHARC : Un ensemble d' outils est nécessaire pour développer les programmes pour la carte SRT-SHARC. Cet ensemble comprend les éléments suivants:
  •
  L'assembleur: Ce logiciel permet l'assemblage des programmes d'utilisation écrit en langage «assembleur ». Il fournit à la sortie un fichier objet qui contient le code bas niveau (code machine) du processeur ADSP21 061-SHARC.
  •
  Le compilateur C : Ce logiciel permet de traduire un fichier source écrit en «C» en langage assembleur. Le fichier ainsi généré sera ensuite assemblé par L'assembleur.
  •
  L'éditeur de lien: Un logiciel qui génère un fichier exécutable par le processeur à partir du fichier déjà assemblé par L'assembleur, et en tenant compte des informations issues du fichier d'architecture.
  •
  SRT_GO: Ce logiciel permet de charger le code bas niveau dans la mémoire interne du processeur et de déclencher l'exécution.
  •
  Le simulateur: il permet de tester sur PC le fonctionnement du programme avant le chargement. Les programmes de la carte SRT-SHARC peuvent être développer en langage assembleur
  ou en C.
  131  Bibliographie
  132 [Asfour-JMR,2000] Asfour. A, Raoof. K, et Fournier. J-M, Nonlinear identification ofNMR spin systems by adaptive filtering, Journal of Magnetic. Resonance, Vol. 145,37-51 (2000).
  [Asfour, 2000]
  Asfour. A, Raoof. K et Fournier. J-M, Identification of NMR spin
  systems by Volterra series, IEEE World Congress on Medical Physics and Biomedical Engineering. Chicago, USA, 23-28 Juillet (2000).
  [Asfour, 2001]
  Asfour. A, Raoof. K, et Fournier. J. M, Filtrage adaptatif RU non
  linéaire: application à l'identification de systèmes à résonance magnétique nucléaire (RMN) , l8/1l Symposium GRETSI'Ol on Signal and Image Processing. Toulouse, France, Septembre (2001).
  [Bellanger, 1989]
  Bellanger. M, Analyse des Signaux et Filtrage Numérique Adaptatif,
  Masson, Paris (1989).
  [Biglierri, 1982]
  E. Biglieri, Theory of Volterra processors and sorne applications,
  Proceedings IEEE ofICASSP. Vol. 82,294-297, (1982).
  [Bittoun, 1984]
  Bittoun. J, Taquin. J, et Sauzade. M, A Computer algorithm for the
  simulation of any nUclear magnetic resonance (NMR) imaging method, Magnetic Resonance in Medicine. Vol. 2, 113-120 (1984).
  [Blütnich, 1982]
  Blütnich. B, et Ziessow. D, Saturation in Hadamard NMR spectroscopy
  and its description by a correlation expansion, Journal of Magnetic Resonance. Vol. 46, 385 -405 (1982).
  [Brenner, 1997]
  Brenner. A .R, Kürsch. J, et Noll. T. G, Distributed large-scale
  simulation ofmagnetic resonance imaging, MAGMA. Vol. 5, 129-138 (1997).
  [Brent, 1986]
  Brent. K et Stewart. B .S, Implementation of a magnetic resonance
  imaging computer simulator and preliminary results, SPIE Medicine. Vol. 626, 200-206 (1986).
  133 [Chester, 1999]
  Chester. D. B, Digital IF filter technology for 3G systems: An
  Introduction, IEEE Communication Magazine, 102-107; Février (1999).
  [Choquet, 1996]
  Choquet. P, Applications Vétérinaires de l'IRM Dédiée à Bas Champ,
  Thèse de Doctorat, Discipline: Imagerie Médicale, Tours, 1996, N°-96TOURS305 .
  [Constantinesco, 1997]
  Constantinesco. A, Choquet. P, et al., Low-Field dedicated and
  desktop magnetic resonance Imaging systems for agricultural and food applications, Magnetic
  Resonance in Chemistry. Vol. 35, 569-575 (1997).
  [Ernst, 1989] Ernst. R. R, Bodenhausen. G, et Wokaun. A, Princip les of Nuc1ear Magnetic Resonance in One and Two Dimensions, Clarendon Press, Oxford (1989).
  [Fakhouri, 1980]
  Fakhouri. S. Y, Identification of the Volterra kernels of nonlinear
  systems, IEE Proceeding, Vol. 127,296-304, (1980).
  Gries. P et Constantinesco. A., Développement d'un imageur de
  [Gries, 1987]
  laboratoire à 0.094T à partir d'un mini spectromètre RMN, Journal de Biophysique et
  Biomécanique. Vol. Il : 139-143 (1987).
  Guinet. CL et Grellet. J, Introduction à l'IRM: De la Théorie à la
  [Guinet, 1992]
  Pratique, Masson, Paris (1992).
  Hoult. D. l, The solution of the Bloch equations in the presence of a
  [Hoult, 1979]
  varying BI field-an approach to selective pulse analysis, Jouranl of Magnetic. Resonnace. Vol. 35 B, 69-86 (1979).
  [Kaulisch, 1996]
  Kaulisch. T, Kolbe. H, et Ziessow. D, Non-linear response theory in nD
  NMR with complex stochastic excitation, Journal of Magnetic Resonance. Series A 121, 4249 (1996).
  [Koh, 1985]
  Koh. T, Powers. E. J, Second-order Volterra filtering and its application to
  nonlinear system identification, IEEE Transaction ASSP, Vol. 33, 1445-1 455, (1985).
  134 [Lawrence, 1981]
  Lawrence. P. J, Estimation of the Volterra functional series of a
  nonlinear system using frequency-response data, IEE Proceeding, Vol. 128,206-210, (1981).
  [Liang, 1999] Liang. Z. P, et Lauterbur. P. C, Princip les of Magnetic Resonance Imaging : a Signal Processing Perspective, IEEE Press, New Press (1999).
  [Mansfield]
  Mansfield. P et Morris. A. P.G, NMR Imaging in Biomedicine, Academic
  Press, New York (1989).
  [Mathews 1991]
  Mathews. V. J, Adaptive polynomial filters, IEEE SP Magazine, Vol.8,
  10-26, July (1991).
  [Michaut, 1992]
  Michaut. M, Méthode Adaptatives pour le Signal: outils mathématiques
  et mise en oeuvre des algorithmes, Hermès, Paris (1992).
  [Morgan 1999]
  Morgan. P. N, Iannuzzelli. R. J, Epstein. F. H, et al., Real-time cardiac
  MRI using DSP's, IEEE transaction on medical Imaging,. Vol. 18,649-653 (1999).
  [Nguyen, 1996]
  Nguyen. V, Système de Création de Champ Magnétique Homogène à
  aimants ferrites pour l'imagerie RMN, Thèse de Doctorat, Discipline: Génie Electrique, Grenoble, 1996.
  [Olsson, 1995]
  Olsson. M. B. E, Wirestam. R, et Presson. B. R, A computer simulation
  for MR imaging: application to RF and static magnetic field imperfections, Magnetic
  Resonance in Medicine, Vol. 34,612-617 (1995).
  [Pimmel, 1990]
  Pimmel. P, Les Antennes en Résonance Magnétique Nucléaire:
  Fonctionnement et Réalisation - Résonateurs pour l'Imagerie et pour la Spectroscopie In
  Vivo, Thèse de Doctorat, Lyon, 1990, N°-101.90.
  [Raoof, 2002]
  Raoof. K, Asfour. A, et Fournier. J. M, A Complete digital MRI system
  at low magnetic field (0.1 Tesla), IEEE Instrumentation and Measurement Technology
  Conference IMTC2002, Anchorage, USA, 20-23 Mai 2002.
  135 [Raoof, 1994]
  Raoof. K et Piatszek. E, A novel real time adaptive filtering system for
  biomedical signal processing, International Conference on Signal Processing and
  Technology, Dallas, USA, Octobre (1994).
  Ravier. S, Gangi. A, Choquet. P et al., IRM Interventionnelle: Intérêt
  [Ravier, 1996]
  du bas champ magnétique et des séquences ultra rapides d'équilibre dynamique 3D pour le suivi d'une aiguille non ferro-magnétique, Journal de Radiologie. Vol. 77 : 105-110 (1997).
  [Ravier, 1994]
  Ravier. S, Développements de Séquences d' Imagerie pour un Prototype
  d' IRM Dédiée à 0.1 Tesla, Thèse de Doctorat, Discipline : Imagerie Médicale, Strasbourg, 1994, N°-1883.
  [Shynk, 1989] Shynk. J. J, Adaptive IIR filtering, IEEE ASSP Magazine, Vo1.6, 4-21, Avril (1989).
  [Stormont]
  Storrnont. R. S et al. , RF synthesizer for an NMR instrument, Brevet Etats
  Unis, United States Patent, No. 4952877, (1990).
  [Summer, 1986]
  Summer. R. M, Axel. L, Israel. S, A computer simulation for nuclear
  magnetic resonance imaging, Magnetic Resonance in Medicine, Vol. 3, 363-376(1986).
  [Walden, 1999]
  Walden. H. R, Performance Trends for Analog-to-Digital Converters, ,
  IEEE Communication Magazine, 96-101, February (1999).
  [Widrow, 1975]
  Widrow. B, et al., Adaptive noise cancelling: princip les and application,
  Proceedings of the IEEE, Vo1.63, No.12, 1692-1716, (1975).
  136 
meta: {'path': '38/tel.archives-ouvertes.fr-tel-00473212-document.txt'}