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

Procede et dispositif d'analyse spectrale en temps reel de signaux instationnaires complexes

Info

Publication number
EP0464115A1
EP0464115A1 EP90905549A EP90905549A EP0464115A1 EP 0464115 A1 EP0464115 A1 EP 0464115A1 EP 90905549 A EP90905549 A EP 90905549A EP 90905549 A EP90905549 A EP 90905549A EP 0464115 A1 EP0464115 A1 EP 0464115A1
Authority
EP
European Patent Office
Prior art keywords
signal
model
parameters
signals
spectral analysis
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
EP90905549A
Other languages
German (de)
English (en)
Inventor
Guy 3 résidence de Chevreuse DEMOMENT
Alain Herment
Claude Arcile
Indira Mouttapa
Amrane Univ.Des Sciences Et De La Houacine
Pierre Peronneau
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institut National de la Sante et de la Recherche Medicale INSERM
Original Assignee
Institut National de la Sante et de la Recherche Medicale INSERM
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institut National de la Sante et de la Recherche Medicale INSERM filed Critical Institut National de la Sante et de la Recherche Medicale INSERM
Publication of EP0464115A1 publication Critical patent/EP0464115A1/fr
Ceased legal-status Critical Current

Links

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/06—Measuring blood flow
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01F—MEASURING VOLUME, VOLUME FLOW, MASS FLOW OR LIQUID LEVEL; METERING BY VOLUME
    • G01F1/00—Measuring the volume flow or mass flow of fluid or fluent solid material wherein the fluid passes through a meter in a continuous flow
    • G01F1/66—Measuring the volume flow or mass flow of fluid or fluent solid material wherein the fluid passes through a meter in a continuous flow by measuring frequency, phase shift or propagation time of electromagnetic or other waves, e.g. using ultrasonic flowmeters
    • G01F1/663—Measuring the volume flow or mass flow of fluid or fluent solid material wherein the fluid passes through a meter in a continuous flow by measuring frequency, phase shift or propagation time of electromagnetic or other waves, e.g. using ultrasonic flowmeters by measuring Doppler frequency shift
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P5/00—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft
    • G01P5/24—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft by measuring the direct influence of the streaming fluid on the properties of a detecting acoustical wave
    • G01P5/241—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft by measuring the direct influence of the streaming fluid on the properties of a detecting acoustical wave by using reflection of acoustical waves, i.e. Doppler-effect
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/02—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems using reflection of acoustic waves
    • G01S15/50—Systems of measurement, based on relative movement of the target
    • G01S15/58—Velocity or trajectory determination systems; Sense-of-movement determination systems
    • G01S15/582—Velocity or trajectory determination systems; Sense-of-movement determination systems using transmission of interrupted pulse-modulated waves and based upon the Doppler effect resulting from movement of targets

Definitions

  • the invention relates to adaptive spectral analysis in real time of unsteady signals with complex values representing a physical phenomenon.
  • the invention can be used whenever it is necessary to analyze in real time such signals representable by a self-regressive parametric model, making the assumption that the coefficients of the model remain stationary over short time intervals. It applies in particular to all fields of medical imaging, for example by NMR, and more generally each time that a spectral analysis must be carried out on organic signals, for example in electrocardiography and electroencephalography. However, it finds a particularly important application in the spectral analysis of the signals provided by ultrasonic Doppler velocimeters used to study the blood flow in the cardiac chambers or the vessels.
  • parametric spectral analysis is relatively little used in Doppler velocimetry and medical ultrasound ultrasound. It remains at the laboratory stage. It consists in modeling the signals using a parametric model, generally of the self-regressive type with adjusted mean, whose parameters are adjusted by optimization using algorithms of the least squares type.
  • This approach has the advantage of making it possible to work on short sections of the signal but it requires the use of low order models and consequently has a limited frequency resolution, which is particularly troublesome for spectral analysis of signals. strongly unsteady.
  • the known methods do not treat demodulated signals with complex values as such, but deal separately with each of the components, respectively in phase and in quadrature.
  • the invention aims in particular to provide a method and a device for spectral analysis which respond better than those previously known to the requirements of practice, in particular in that they improve the frequency resolution, allow a more precise analysis of the signals and restrict the variability of the spectra obtained to represent flows under similar conditions.
  • J 1 (a) is the sum of the least squares [y (n) -z (n)] 2 , z denoting the sampled function modeling the signal via the parameters a to be determined and y the function represented by the samples of the signal and where ⁇ (a) is a regularizing functional taking into account a priori knowledge of the nature of the function or obtained by optimization.
  • the method makes it possible to directly obtain parameters for calculating the power spectral density of the signal in the analysis window.
  • the invention also provides a spectral analysis device directly adaptable to a conventional Doppler velocimeter, without modification of the existing probe and its associated electronics, or of the display means.
  • FIG. 1 is a block diagram showing the parameters involved in Doppler velocimetry of blood flow in the particular case of a continuous Doppler technique
  • FIG. 2 is a diagram showing the envelope of a representative Doppler spectrum
  • FIG. 3 is a block diagram of the demodulation of the Doppler signal
  • FIG. 4 and 5 are block diagrams showing a possible implementation of the method according to the invention.
  • the non-stationarity of the signal is sufficiently slow compared to the speed of variation of the amplitude of the signal, the latter can be considered as "locally stationary", that is to say stationary on a slice of analysis of N samples, N typically being of the order of a few tens.
  • y p (n) t [y (np), ..., y (n-1)] where a t is the vector of the coefficients of the model whose order is p, and where b (n) represents the process signal generator, with variance ⁇ 2 unknown.
  • the spectral analysis is then carried out according to a criterion which is often that of the least squares, consisting in minimizing J 1 (a): J 1 (a) - z (n)] 2 (4)
  • J 2 (a) J 2 (a) + ⁇ (a) (5)
  • ⁇ (a) is a regularizing functional, of which there are many examples.
  • ⁇ is therefore a scalar, product of the transposed matrix
  • the transpose is obviously a line vector.
  • the choice of the regularizing functional ⁇ (a) translates the a priori information of the nature of the local spectral properties of the signal analyzed. This information is that the spectral density of unknown power is a "soft" function, that is to say having a certain regularity, compared to an unorganized spectrum made up of a multitude of lines without relation between them. For the sake of convenience implementation, a quadratic functional is chosen, which has the advantage of leading to an explicit and linear solution with respect to the observed values y (n).
  • H (f) is the frequency transfer function of the whitening filter of the analyzed signal, which is also the inverse of its generator filter which is described by equations (1) and (2).
  • D K signifies a non-soft power spectral density, all the more so as the order k of the derivative is high.
  • the invention uses an approach which can be considered to be in two stages, using in particular the observations that we very often deal in practice with signals for which we can carry out a local study of the phenomenon by dividing the time horizon into blocks of length N, and that we can assume a local stationarity in each block, so as to make it possible to exploit local invariance properties of the model, by therefore carrying out an adaptive spectral analysis reflecting a slow variation of the parameters a -vis signal amplitude fluctuations, namely:
  • This model is sought for the value of the regularization coefficient ⁇ and that of the variance ⁇ 2 of the generating process which maximizes the likelihood of the model.
  • the first step involves organizing the samples of the signal into consecutive blocks of length N and performing regularization on each block.
  • Y m [y (-p + 1), y (-p + 2), ..., y (0), y (1), ..., y (mp)] t
  • a condensed writing of the algorithm is obtained by projection, by only showing the components actually intervening at each recurrence in the calculation.
  • k m (i) [O t i + 1 , k (i) t , O t mpi-1 ] t
  • This algorithm has the advantage of applying both to "pre-priestly” problems (treatment of the stationary case or of an isolated block in the non-stationary case) and to problems of the covariance type (treatment of adjacent blocks in the non-stationary case).
  • ⁇ P m (2) DP m (1/0) D t -k m (1) r (1) -1 k m ( 1) * - P m (1/0) with
  • the initial factorization not being unique, one chooses advantageously for M (1) the signature matrix of ⁇ P m (2) which provides the beginning of an algorithm in square root whose numerical qualities are better.
  • ⁇ P m (2) - ⁇ O S m S * m + ⁇ O v m v * m - k m (1) r (1) -1 k m (1) *, of rank at most equal to 3 since it breaks down into a sum of three dyadic or anti-scalar matrices, each of rank equal to 1.
  • ⁇ Pm (2) in the form:
  • the noise generating the observed signal b (n) which is also complex, has real and imaginary parts whose a priori probability laws are also normal, centered, independent, and of variance (1/2) ⁇ 2 .
  • pre-winding A solution to simplify the algorithm can be described as "pre-winding". It assumes that the p values prior to the initial instant in the measurement window are zero:
  • ⁇ P m (2) - ⁇ O s m s * m + ⁇ O v m v * m
  • the complexity of the calculation is thus significantly reduced.
  • the hyperparameters can be chosen a priori. Otherwise, a calculation of the likelihood of the hyperparameters is useful.
  • Kalman filter makes it possible to calculate by recurrence the likelihood of the hyperparameters (or its logarithm) using the quantities involved in the algorithm.
  • the model to be retained is the one which maximizes this likelihood compared to hyperparameters.
  • the method is adaptive since the choice of a priori values may depend on the data.
  • V ( ⁇ , ⁇ 2 / y) f [y (1)] f [y (n) / y (1), ..., y (n-1)]
  • the complete process includes, when determining the optimal hyperparameters:
  • the speed field studied is illuminated by a monochromatic or pulsed acoustic wave coming from the transmitter 10 and the wave reflected by the moving material targets (red cells) is measured. under the action of this speed field, collected by a receiver 12 comprising a demodulator.
  • the received signal r (t) would be deduced from the transmitted signal a (t) by a simple frequency translation. It is not the case, because of the contribution of a large number of elementary targets having a dispersion compared to an average ideal target, the received signal has a continuous spectrum, of narrow bandwidth compared to the frequency of the carrier (a few kHz compared to a few MHz). In the frequent case of a Doppler signal with pulsed emission, the shape of the emitted wave also comes into play. In both cases, a narrow band spectrum, of the kind shown in FIG. 2, is obtained, limited to two symmetrical bands of width 2F centered at -f 0 and f 0 respectively, F being much lower than the frequency f 0 of the carrier.
  • demodulation generally synchronous, which leads to decomposing the received signal r (t) into two components, respectively in phase p (t) and in quadrature q (t).
  • These low frequency signals (a few kHz) require only a reduced sampling rate because their useful spectrum is reduced by a factor of approximately 1000 relative to the carrier.
  • the assembly can be that shown diagrammatically in FIG. 3.
  • p (t) and q (t) are low frequency signals generated by the demodulator 14 of FIG. 3, constituting respectively the real part and the imaginary part of the "Doppler signal” with complex values y (t).
  • the invention therefore deals simultaneously with p (t) and q (t), which is possible because the spectral analysis method described above was precisely designed to apply to signals with complex values.
  • the initialization mode must also be chosen. It is defined by:
  • a 0 can either be zero, or the result of the processing of the previous block; the choice between the version known as "covariance” and that known as “prefenêtrée”.
  • N, p and ⁇ results from a compromise. Increasing N and p improves the quality of the model, but we are limited in this way by the non-stationarity of the signal which pushes to reduce N and by the computational load which pushes to reduce p.
  • the initialization can be done by pre-windowing or by covariance, as indicated above.
  • the pre-windowing currently used, decreases the number of multiplications to be performed, but harms the likelihood of the model, especially if p is close to N.
  • the choice of the hyperparameter ⁇ is essential because it regulates the degree of softness of the solution and must maximize the likelihood of the signal model compared to the data: the values of the calculated AR parameters do not in fact depend on the variances ⁇ 2 and t but only from their ⁇ ratio.
  • the likelihood variation curve as a function of ⁇ varies very little in the vicinity of the maximum and that the results of the analysis are little affected by variations of an order of magnitude of this parameter.
  • the variations in the optimal value of ⁇ from one data block to another are also of an order of magnitude.
  • the complexity of the method can therefore be reduced by substituting, for the maximization stage of the likelihood of the model which must be performed on each data block, a prior step of classification of the signals according to the constant values to be assigned to ⁇ . For this to be possible, it is essential to normalize ⁇ with respect to the strength of the Doppler signal. Indeed, the coefficient intervenes in the minimized criterion:
  • J 1 (a) has the dimension of the energy of the signal and the parameters AR are dimensionless.
  • the coefficient ⁇ therefore has a dimension and varies like the square of the scale in which the Doppler signal is measured.
  • ⁇ ⁇ 0 ⁇ 2 (N / p)
  • ⁇ 0 the new adjustment parameter, this time dimensionless
  • ⁇ 2 the variance of the signal which can be supplied by the acquisition module of the device used (24).
  • the invention is susceptible of numerous variant embodiments.
  • the device used for the implementation of the The process may have the general constitution shown in FIG. 4, which seeks a compromise between the quality of the results and the volume of the calculations, is adaptable to existing devices and uses the intrinsic possibilities of digital signal processing processors currently available.
  • the device has two processors coupled to perform the Doppler spectrum calculations.
  • the first, 20, calculates the parameters of a model of the received Doppler signal, using the B-CAR algorithm.
  • the second, 22, provides the spectrum from this model.
  • This second processor 22 controls an acquisition card 24 which acquires the two Doppler channels simultaneously, at a sampling frequency which is for example 40 kHz for the study of blood circulation, chosen to optimize the compromise between the duration of the slices analyzed and the frequency accuracy of the spectral analysis.
  • a single converter is used.
  • the card can be provided to detect any saturation of the Doppler signals.
  • the processor 20 is programmed to scale the hyperparameter ⁇ of the method from the estimation of the instantaneous power of the Doppler signal supplied by the processor 22, and calculate the AR parameters of the Doppler signal on each data block transmitted by the same module.
  • the processor 22 calculates the spectrum of the signal on a block, from the parameters AR supplied by the processor 20 under the control of the display module 26 which fixes the number of points on the frequency axis and the number of gray levels; provision may be made for calculating a possible continuous component in the acquisition chain. It must also evaluate the instantaneous signal power, used by the processor 20 to optimize the results provided.
  • the architecture of the device can be that shown in FIG. 5.
  • the processor 22 executes the algorithm for calculating the power spectral density and manages the data input and output interfaces 28-30-32.
  • the processor 20 executes the B-CAR algorithm.
  • a mini host computer constitutes a master machine which loads the program into the processors and which controls the execution of the programs.
  • a ROM implementation of the program would avoid the need for a mini host computer.
  • the acquisition 30 and display 32 interfaces transfer the data, on the one hand between the acquisition device and the processor 22, on the other hand between the display device 26 and this same processor 22.
  • Each processor is provided with 'own working memory 36.

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Multimedia (AREA)
  • Biomedical Technology (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Hematology (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Biophysics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Electromagnetism (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Complex Calculations (AREA)

Abstract

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

Description

Procédé et dispositif d'analyse spectrale en temps réel de signaux instationnalres complexes
L'invention concerne l'analyse spectrale adaptative en temps réel de signaux instationnalres à valeurs complexes représentant un phénomène physique.
L'invention est utilisable chaque fois qu'il est nécessaire d'analyser en temps réel de tels signaux représentables par un modèle paramétrique auto-régressif, en faisant l'hypothèse que les coefficients du modèle restent stationnaires sur des intervalles de temps courts. Elle s'applique notamment à tous les domaines de l'imagerie médicale, par exemple par RMN, et plus généralement chaque fois qu'une analyse spectrale doit être effectuée sur des signaux organiques, par exemple en électrocardiographie et électroencéphalographie. Elle trouve toutefois une application particulièrement importante dans l'analyse spectrale des signaux fournis par les vélocimètres Doppler ultra-sonores servant à étudier l'écoulement sanguin dans les cavités cardiaques ou les vaisseaux.
A l'heure actuelle, l'analyse spectrale paramétrique est relativement peu utilisée en vélocimétrie Doppler et en échographie ultra-sonore médicale. Elle reste au stade du laboratoire. Elle consiste à modéliser les signaux à l'aide d'un modèle paramétrique, généralement du type auto-régressif à moyenne ajustée, dont les paramètres sont ajustés par optimisation à l'aide d'algorithmes du type moindres carrés. Cette approche a l'avantage de permettre de travailler sur des tranches courtes du signal mais elle impose d'utiliser des modèles d'ordre faible et a en conséquence une résolution en fréquence limitée, ce qui est particulièrement gênant pour l'analyse spectrale de signaux fortement instationnalres. Par ailleurs, les méthodes connues ne traitent pas les signaux démodulés à valeurs complexes comme tels, mais traitent séparément chacune des composantes, respectivement en phase et en quadrature.
Les appareils du commerce utilisent des méthodes d'analyse spectrale non-paramétriques du type analyse de Fourier pour pouvoir en effectuer la mise en oeuvre à l'aide d'algorithmes de transformation de Fourier rapide. Mais ces méthodes nécessitent de travailler sur des tranches longues du signal pour obtenir une estimation convenable de la représentation spectrale recherchée, ce qui est incompatible avec la nature des instationnarités rencontrées dans les domaines précédemment cités.
L'invention vise notamment à fournir un procédé et un dispositif d'analyse spectrale répondant mieux que ceux antérieurement connus aux exigences de la pratique, notamment en ce qu'ils améliorent la résolution en fréquence, permettent une analyse plus précise des signaux et restreignent la variabilité des spectres obtenus pour représenter des écoulements dans des conditions voisines.
Dans ce but, l'invention propose notamment un procédé d'analyse spectrale adaptative par modélisation du phénomène représenté par le signal, caractérisé en ce qu'on modèlise des tranches successives du signal échantillonné, par un processus auto-régressif d'ordre élevé, proche du ou égal au maximum autorisé par la durée des tranches et la fréquence d'échantillonnage, en faisant une estimation des paramètres caractérisant le modèle par filtrage de Kalman rapide utilisant un algorithme qui est simplifié dans toute la limite autorisée par la stationnarité locale du modèle dans la tranche et en ce qu'on détermine la densité spectrale de puissance à partir des paramètres estimés ; pour éviter les risques d'instabilité liés au choix d'un ordre élevé, on régularise l'estimation des paramètres du modèle auto-régressif par mise en oeuvre d'une contrainte consistant à minimiser un critère de la forme : J2(a) = J1(a) + Ω(a)
où J1(a) est la somme des moindres carrés [y(n)-z(n)]2, z désignant la fonction échantillonnée modélisant le signal par l'intermédiaire des paramètres a à déterminer et y la fonction représentée par les échantillons du signal et où Ω(a) est une fonctionnelle régularisante tenant compte de la connaissance a priori de la nature de la fonction ou obtenue par optimisation.
Grâce à l'utilisation d'un modèle d'ordre élevé, on obtient une très bonne définition en fréquence. Le procédé permet d'obtenir directement des paramètres de calcul de la densité spectrale de puissance du signal dans la fenêtre d'analyse. On peut notamment utiliser un ordre de modèle égal au nombre de points d'échantillon- nage dans le signal Doppler, ce qui donne une excellente définition en fréquence, alors que les méthodes actuelles d'analyse spectrale par modèle auto-régressif et critère des moindres carrés ont un ordre de grandeur 3 à 7 fois inférieur au nombre de points échantillonnés.
L'invention propose également un dispositif d'analyse spectrale directement adaptable sur un vélocimètre Doppler classique, sans modification de la sonde existante et de son électronique associée, ni des moyens de visualisation.
L'invention sera mieux comprise à la lecture de la description qui suit du procédé et du dispositif selon l'invention et du cheminement qui a permis d'y aboutir. Elle se réfère aux dessins qui l'accompagnent, dans lesquels :
- la Figure 1 est un schéma de principe montrant les paramètres qui interviennent en vélocimétrie Doppler du flux sanguin dans le cas particulier d'une technique de Doppler continu ;
- la Figure 2 est un schéma montrant l'enveloppe d'un spectre Doppler représentatif ;
- la Figure 3 est un schéma de principe de la démodulation du signal Doppler ; - les Figures 4 et 5 sont des synoptiques montrant une implémentation possible du procédé suivant l'invention.
Modélisation des signaux non-stationnaires
Avant d'aborder la description de l'invention, il convient de rappeler des indications sur la modélisation des signaux représentant des phénomènes non stationnaires et sur la régularisation du modèle.
Dans le cas où la non-stationnarité du signal est suffisamment lente devant la vitesse de variation de l'amplitude du signal, celui-ci peut être considéré comme "localement stationnaire", c'est-à-dire stationnaire sur une tranche d'analyse de N échantillons, N étant typiquement de l'ordre de quelques dizaines. Les échantillons observés y(n), n = 1, 2,..., N, sont décrits par un modèle auto-régressif (AR) exprimé par : y(n) = z(n)+b(n) (1)
z(n) = (n-i) = atyp(n) (2)
avec at - [ap,ap-1,..., a1] (3) yp(n)t = [y(n-p),..., y(n-1)] où at est le vecteur des coefficients du modèle dont l'ordre est p, et où b(n) représente le processus générateur du signal, de variance σ2 inconnue.
L'analyse spectrale s'effectue alors suivant un critère qui est souvent celui des moindres carrés, consistant à minimiser J1(a) : J1(a) - z(n)]2 (4)
Le modèle étant linéaire, les algorithmes de calcul sont rapides. Mais l'obtention d'une bonne approximation est subordonnée au choix d'un ordre p aussi élevé que possible ; or, si on prend p proche de N ou égal à lui, la minimisation de J1 n'est pas un critère satisfaisant par suite de la variance excessive de la solution.
Pour écarter cette difficulté et stabiliser la solution, on a proposé de régulariser le modèle en ajoutant à J1 une fonctionnelle stabilisatrice, ce qui permet de conserver une valeur élevée de p (proche de ou même égale à N). Le critère consiste alors à minimiser J2(a) :
J2(a) = J2(a) + Ω(a) (5) où Ω(a) est une fonctionnelle régularisante, dont il existe de nombreux exemples.
Cette fonctionnelle Ω sera cependant en règle générale de la forme :
Ω(a) = |a-a0|t M (a-a0)
Ω est donc un scalaire, produit de la matrice transposée |a-a0|t, d'une matrice symétrique et définie positive de dimension p p (p étant la dimension de a) et (a-a0). La transposée est évidemment un vecteur ligne.
Si M est la matrice identité, on arrive à la forme simplifiée qui sera donnée en (7) plus loin.
Le choix de la fonctionnelle régularisante Ω(a) traduit l'information a priori de la nature des propriétés spectrales locales du signal analysé. Cette information est que la densité spectrale de puissance inconnue est une fonction "douce", c'est-à-dire présentant une certaine régularité, comparée à un spectre inorganisé constitué d'une multitude de raies sans relation entre elles. Pour des raisons de commodité de mise en oeuvre, on choisit une fonctionnelle quadratique, ce qui présente l'avantage de conduire à une solution explicite et linéaire par rapport aux valeurs observées y(n).
Une mesure de "douceur spectrale" est constituée par Dk = df k = 0,1,2,...
où H(f) est la fonction de transfert en fréquence du filtre blanchisseur du signal analysé, qui est aussi l'inverse de son filtre générateur qui est décrit par les équations (1) et (2). Une valeur élevée de DK signifie une densité spectrale de puissance non douce, et ce d'autant plus que l'ordre k de la dérivée est élevé. On utilise une douceur spectrale d'ordre zéro
DQ = | H(f) |2 df : 1 + a2 m
qui est, à une constante près, la norme euclidienne du vecteur des paramètres a, ce qui conduit dans la plupart des cas à adopter la fonctionnelle régularisante
Ω(a) = m (6)
où le coefficient μ est un nombre réel positif appelé "coefficient de régularisation" et dont il sera vu qu'il joue un très grand rôle dans le fonctionnement de l'invention et qui est choisi par une méthode de maximum de vraisemblance a posteriori.
Modélisation à régularisation suivant l'invention.
L'invention utilise une démarche qui peut être considérée comme en deux étapes, utilisant notamment les constatations que l'on a très souvent affaire, dans la pratique, à des signaux pour lesquels on peut effectuer une étude locale du phénomène par division de l'horizon temporel en blocs de longueur N, et qu'on peut faire l'hypothèse d'une stationnarité locale dans chaque bloc, de manière à permettre d'exploiter des propriétés d'invariance locale du modèle, en effectuant donc une analyse spectrale adaptative reflétant une variation lente des paramètres a vis-à-vis des fluctuations d'amplitude du signal, à savoir :
1. On recherche un modèle localement stationnaire, c'est-à-dire un modèle invariant équivalent défini par l'optimum du critère (5) calculé sur une fenêtre temporelle réduite, en effectuant un compromis entre une trop grande longueur de fenêtre (qui ne permet pas de suivre fidèlement les non-stationnarités du signal) et une trop courte longueur (qui ne permet pas d'estimer le modèle dans de bonnes conditions), ce qui conduit à choisir un modèle AR de l'ordre p le plus élevé possible, proche de N ou égal à N, imposant de stabiliser le modèle.
Ce modèle est recherché pour la valeur du coefficient de régularisation μ et celle de la variance σ2 du processus générateur qui rendent maximale la vraisemblance du modèle.
2. On simplifie les calculs de maximisation de vraisemblance du modèle en adoptant une solution sous-optimale :
- on réduit le domaine de variation de μ de l'équation (6) du critère de minimisation à un ensemble fini et présélectionné de valeurs discrètes, et
on calcule, par un algorithme qui peut être rapide, la solution optimale à et éventuellement la vraisemblance correspondante V(σ2, μ/y) à partir des échantillons du signal observé, pour retenir la plus élevée.
La première étape conduit à organiser les échantillons du signal en blocs consécutifs de longueur N et à effectuer sur chaque bloc une régularisation. On impose une contrainte liant localement l'ensemble des paramètres et traduisant une information de douceur spectrale : cela conduit à minimiser, sur chacun des blocs, un critère de la forme :
J2(a)= J1(a) + μ |a-a0 |2 (7)
Lorsque le processus générateur du signal analysé est supposé gaussien, cela revient à admettre que les paramètres a sont caractérisés, dans une interprétation bayésienne, par la loi de probabilité a priori :
L'initialisation du vecteur des paramètres AR se fait en choisissant a0 = 0 dans le cas d'un signal stationnaire. On verra plus loin que cette façon de procéder peut être conservée pour certains signaux non stationnaires. Dans les autres cas on choisira, pour chaque bloc courant, le résultat du traitement du bloc précédent lorsque les blocs sont adjacents.
Il reste encore à fixer la valeur de μ et de σ 2 pour pouvoir effectuer la minimisation du critère (7). Ce choix est fait, comme il a été dit plus haut, en maximisant la vraisemblance du modèle. On verra que μ etσ2, appelés hyperparamètres, peuvent être découplés dans cette maximisation.
Pour exposer la première étape, il convient de donner les équations d'état associées au problème et de montrer comment on peut substituer au filtre de Kalman standard (qui fait appel à une équation de Riccati n'exploitant pas la propriété de stationnarité locale du modèle qui es base des méthodes adaptatives, ni la propriété de age des vecteurs yp(n) successifs qui est la conséquence du choix d'un modèle AR) un filtre plus simple, mettant en oeuvre les équations de Chandrasekhar, pour le calcul des paramètres â qui minimisent, μ et σ2 étant fixés, le critère (7), ainsi que pour le calcul de la vraisemblance correspondante. On verra que ce filtre, à modèle d'état invariant, permet d'introduire la contrainte de régularisation.
Pour cela, on fait apparaître explicitement ces propriétés en définissant un "vecteur-paramètre étendu" ou vecteur d'état à m composantes : am(i) - [Oti-1, ap,ap-1,..., a1,0,0,..., 0]t, avec m>p et on considère le vecteur global des données :
Ym = [y(-p+1),y(-p+2),..., y(0), y(1),..., y(m-p)]t
Le modèle linéaire défini par les équations (1) et (2) peut alors se mettre sous la forme du modèle d'état équivalent suivant : am(i+1) = Dam(i) i = 1,2,... y(i) = γtm am(i)+b(i) où D est l'opérateur de décalage
Le calcul récurrent de la solution peut alors être effectué par le filtre de Kalman suivant : âm(i+1/i)=Dâm(i/i-1) +km(i)r(i)-1[y(i)-yt mam(i/i-1)] (8)
avec les conditions initiales :
où désigne le vecteur simplement conjugué du vecteur ym et k* désigne la matrice transposée conjuguée de k.
En apparence, l'introduction de ces variables supplémentaires est défavorable, car elle augmente les dimensions du vecteur d'état. Mais cette impression est trompeuse, car cette introduction fait apparaître explicitement les propriétés d'invariance que l'on exploite.
On peut d'ailleurs vérifier qu'à chaque récurrence, le nombre de coordonnées non-nulles du vecteur gain km(i) reste égal à p.
Ce modèle d'état étant invariant (D et ym étant constant) et à variance stationnaire (σ2 étant constant sur un bloc), on obtient directement un algorithme rapide en remplaçant les équations standard par des équations de Chandrasekhar. Non seulement ces équations permettent une résolution rapide, mais encore elles permettent de passer, de façon naturelle, à des algorithmes en racine carrée offrant une bonne stabilité numérique.
On désignera par la suite par "Algorithme B-CAR" les règles de régularisation mettant en oeuvre l'algorithme, appliquées à des blocs de taille appropriée à la nature du signal. Ces équations ne propagent que les incréments des quantités nominales de l'équation (9), qui sont eux-mêmes factorisés. L'application directe de ces techniques de factorisation au filtre (8) (9) conduit aux équations suivantes :
On utilise dans ces équations le fait que M et Pm sont des matrices à symétrie hermitienne : M*=M et P*m=Pm. Une écriture condensée de l'algorithme est obtenue par projection, en ne faisant apparaître que les composantes intervenant effectivement à chaque récurrence dans le calcul. On définit pour cela un vecteur-gain effectif k(i), de dimension p, contenant les coordonnées non-nulles de km(i) : km(i) = [Ot i+1, k(i)t, Ot m-p-i-1]t
De la même façon, si α est la dimension de M(i), on peut aussi écrire :
B(i) = [Oα,k, S(i)t, Oα,m-p-i-1]t où S(i) est une matrice de dimensions [(p+1),α]. L'algorithme peut alors se mettre sous la forme
1
la remise à jour des paramètres se faisant selon : â(i+1/i)-â(i/i-1)+k(i)r(i)-1[y(i)-yp(i)tâ(i/i-1)] (13)
Cet algorithme a l'avantage de s'appliquer aussi bien à des problèmes "préfenêtrés" (traitement du cas stationnaire ou d'un bloc isolé dans le cas non-station- naire) qu'à des problèmes du type à covariance (traitement de blocs adjacents dans le cas non-stationnaire).
L'initialisation de l'algorithme est conditionnée par la matrice de covariance a priori P (1/0) dont se déduisent B(l) et M(l) lors de la factorisation de l'incrément initial : δPm(2) = PM(2:1)-Pm(1/0) - B(1)M(1)B(1)t
L'incrément initial de covariance peut également s'écrire, en revenant aux équations sous forme étendue : δPm(2) = DPm(1/0) Dt-km(1)r(1)-1km(1)* - Pm(1/0) avec
Pm(1/0)
Dans le cas général, si l'on pose
P(1/0)= v0 alors on peut décomposer DPm( 1/0 ) D t-P m( 1/0 ) = -ν0( sm s* m sm g* m + gm s* m ) + ν0( vm v* m + vm d* m + dm v* m ) où
La matrice P(1/0) est choisie diagonale, et le rang de δmP(2) est en général α=3. La factorisation initiale n'étant pas unique, on choisit avantageusement pour M(1) la matrice de signature de δPm(2) qui fournit l'amorce d'un algorithme en racine carrée dont les qualités numériques sont meilleures.
Si la matrice de covariance a priori P(1/0) est diagonale, les vecteurs gm et dm sont identiquement nuls et l'incrément initial δPm(2) se réduit à :
δPm(2) = -νOSm S* m + νO vm v* m - km(1)r(1)-1km(1)*, de rang au plus égal à 3 puisqu'il se décompose en une somme de trois matrices dyadiques ou antiscalaires, chacune de rang égal à 1. On peut écrire δPm(2) sous la forme :
δPm( 2)= B(1) M(1) B(1) *
avec comme facteur central
M(1) =
qui conduit, une fois la projection effectuée pour obtenir l'algorithme en dimension réduite, à :
sous les conditions, généralement remplies, que :
- les coefficients ai du modèle, qui sont complexes, aient des parties réelles et imaginaires dont les lois de probabilité a priori sont normales, centrées, indépendantes, et de variance (1/2) 2 (on retrouve ainsi la matrice de covariance a priori P( 1/0 )
- le bruit générateur du signal observé b(n), qui est aussi complexe, ait des parties réelle et imaginaire dont les lois de probabilité a priori sont aussi normales, centrées, indépendantes, et de variance (1/2)σ2.
Cette solution rigoureuse utilise les valeurs effectivement observées antérieures à l'instant initial de la fenêtre dans le vecteur général ym. Cette méthode, que l'on peut qualifier de "covariance" suppose que l'on travaille :
- soit sur un bloc de données isolé de longueur double de l'ordre p choisi,
- soit sur un bloc de données adjacent à un bloc antérieur dont le contenu a été conservé.
Le décompte des opérations arithmétiques élémentaires, limité aux seules multiplications, montre que la complexité de cet algorithme est θ(11p) par récurrence. Si on choisit p=N, on voit que la complexité totale pour traiter N données est θ(N2), ce qui est normal pour un algorithme de ce type.
Une solution permettant de simplifier l'algorithme peut être qualifiée de "préfenêtrage". Elle fait l'hypothèse que les p valeurs antérieures à l'instant initial dans la fenêtre de mesure sont nulles :
avec pour conséquences :
δPm(2) = -νOsms*m+νOvmv*m
Le rang de l'incrément initial est α=2. On peut choisir comme facteur initial :
M(1) =
S(1) =
On réduit ainsi sensiblement la complexité du calcul. Lorsque la nature du phénomène représenté par le signal est bien connue, les hyperparamètres peuvent être choisis a priori. Dans le cas contraire, un calcul de la vraisemblance des hyperparamètres est utile.
L'utilisation d'un filtre de Kalman permet de calculer par récurrence la vraisemblance des hyperparamètres (ou son logarithme) à l'aide des quantités intervenant dans l'algorithme.
La recherche du maximum de vraisemblance revient en fait à un problème de moindres-carrés régularisés déjà représentés par l'équation (7). Ce problème a une interprétation bayésienne simple puisque la minimisation de (7) est équivalente à la maximisation du critère de vraisemblance V(a) : V(a) = exp (15)
Lorsqu'on fait intervenir les hyperparamètres μ, σ2, le calcul montre que la solution â est en fait la moyenne de la loi a posteriori définie par la distribution conditionnelle des données :
et la loi a priori :
On peut donc considérer que cette loi a priori spécifie une classe d'estimateurs par l'intermédiaire des paramètres μ et σ2. Ce sont les hyperparamètres du problème. Comme les lois sont normales, la vraisemblance des hyperparamètres est représentée par une intégrale : V(σ2,μ,y) = f(y/a,σ2)f(a/μσ2)da (16)
Le modèle à retenir est celui qui maximise cette vraisemblance par rapport aux hyperparamètres. En ce sens, la méthode est adaptative puisque le choix de valeurs a priori peut dépendre des données.
Enfin, par application séquentielle de la règle de Bayes à la relation, on obtient :
V(μ, σ2/y) = f[y(1)] f[y(n)/y(1),..., y(n-1)]
On utilise la densité marginale conditionnelle de y(n) que permettent d'obtenir les échantillons observés y(1), y(2),...y(n-1) :
Le calcul de cette densité conditionnelle exige de connaître σ2 et a. Mais y(n) = yp(n)tâ(n/n-1)+e(n) où â(n/n-1) est l'estimée optimale (c'est-à-dire conditionnelle au passé (y(1), y(2),..., y(n-1)), et où e(n) est l'innovation du processus observé. On a donc aussi : f[y(n) / y(1),..., y(n-1)]
Reste à évaluer la variance de l'innovation σe(n)2. Mais un filtre de Kalman n'est sensible qu'à la seule variance relative du bruit d'observation (c'est- à-dire ici du processus générateur) et du bruit d'état (c'est-à-dire ici de la variance a priori des paramètres). Le calcul récurrent de â(n/n-1) ne dépendant donc que de μ et pas de la valeur de σ2, on peut prendre arbitrairement σ2 = 1 dans l'algorithme ; on a alors : σe(n)2 = σ2r(n) = σ2 Pm(n/n-1)ym+1]
d'où l'on tire : f[y(n) / y(1),..., y(n-1)]
Cette formule montre que les deux hyperparamètres sont découplés. L'estimée de σ2 selon le maximum de vraisemblance (marginale) est donc :
et la vraisemblance (marginale) de μ est alors : L(μ/y) = + In|r(i)|] + cste
En résumé le procédé complet comprend, lorsqu'on doit déterminer les hyperparamètres optimaux :
- pour chacune des valeurs discrètes de μ choisies pour le calcul de J2, l'estimation d'un vecteur â(N/N) avec l'algorithme (10) (11), avec σ2 = 1 ; le choix des paramètres à et des hyperparamètres σ2 et μ correspondant à la vraisemblance la plus forte ; et le calcul de la densité spectrale de puissance g(f) correspondante à l'aide de la relation classique
g(f) = σ2/ [ ak exp(-2iwfk)|2] (17)
pour
- 1/2 ≤ f ≤ + 1/2
Application de la modélisation suivant l'invention aux signaux de vélocimètrie sanguine Doppler
Le choix de la durée des fenêtres ou blocs de modèle auto-régressif, du nombre d'échantillons à prélever par bloc, de l'ordre du modèle et du critère de régularisation implique de connaître les caractères principaux des signaux Doppler. Ces signaux sont obtenus par un montage du genre schématisé en Figure 1. Le champ de vitesse étudié est éclairé par une onde acoustique monochromatique ou puisée provenant de l'émetteur 10 et on mesure l'onde réfléchie par les cibles matérielles (hématies) se déplaçant sous l'action de ce champ de vitesse, recueillie par un récepteur 12 comportant un démodulateur.
Si le champ de vitesse était stationnaire et uniforme, le signal reçu r(t) se déduirait du signal émis a(t) par une simple translation en fréquence. Ce n'est pas le cas, du fait de la contribution d'un grand nombre de cibles élémentaires présentant une dispersion par rapport à une cible idéale moyenne, le signal reçu a un spectre continu, de largeur de bande étroite par rapport à la fréquence de la porteuse (quelques kHz par rapport à quelques MHz). Dans le cas fréquent d'un signal Doppler à émission puisée, la forme de l'onde émise entre également en jeu. On obtient dans les deux cas un spectre à bande étroite, du genre montré en figure 2, se limitant à deux bandes symétriques de largeur 2F centrées en -f0 et f0 respectivement, F étant très inférieur à la fréquence f0 de la porteuse.
Pour simplifier le traitement des signaux, les appareils existants utilisent une démodulation, généralement synchrone, qui conduit à décomposer le signal reçu r(t) en deux composantes, respectivement en phase p(t) et en quadrature q(t). Ces signaux à basse fréquence (quelques kHz) n'exigent qu'une cadence d'échantillonnage réduite car leur spectre utile est réduit par un facteur 1000 environ par rapport à la porteuse. Le montage peut être celui schématisé en figure 3.
En représentation complexe, on peut écrire :
e( t ) =
où ω0 est la pulsation de la porteuse et où e(t) désigne l'enveloppe complexe de modulation du signal émis, et r(t) = p(t)cos(ω0t)+ q(t)sin(ω0t)
où p(t) et q(t) sont des signaux à basse fréquence générés par le démodulateur 14 de la figure 3, constituant respectivement la partie réelle et la partie imaginaire du "signal Doppler" à valeurs complexes y(t).
Le calcul montre qu'on peut obtenir le spectre du signal reçu r(t) par transformation de Fourier de la covariance du signal Doppler y(t). Par ailleurs p(t) et q(t) sont en général corrélées, ce qui interdit de faire l'hypothèse - pourtant admise dans le traitement habituel du signal Doppler - que la covariance mutuelle RPQ de p et q ou leur densité interspectrale de puissance GPQ(f) est nulle. Or cette hypothèse est nécessaire si on traite p(t) et q(t) séparément.
L'invention traite donc simultanément p(t) et q(t) ce qui est possible du fait que la méthode d'analyse spectrale décrite plus haut a été justement conçue pour s'appliquer à des signaux à valeurs complexes. Choix des paramètres de réglage dans l'application à la vélocimétrie Doppler
Trois paramètres de réglage doivent être choisis, soit a priori (ce qui sera le cas général étant donné que l'on connaît la nature du signal Doppler), soit à l'issue d'une procédure d'optimisation.
Il s'agit de :
- N, durée d'une fenêtre d'analyse ou longueur d'un bloc de données ;
- p, ordre du modèle AR adopté pour décrire le signal Doppler,
- μ, coefficient de régularisation ou hyperparamètre.
Le mode d'initialisation doit aussi être choisi. II est défini par :
le choix des paramètres initiaux a0 au début de chaque bloc: a0 peut être ou bien nul, ou bien le résultat du traitement du bloc précédent ; le choix entre la version dite "covariance" et celle dite "préfenêtrée".
Le choix des trois paramètres N, p et μ résulte d'un compromis. Augmenter N et p améliore la qualité du modèle, mais on est limité dans cette voie par la non- stationnarité du signal qui pousse à réduire N et par la charge de calcul qui pousse à réduire p.
On obtient de bons résultats avec une fréquence d'échantillonnage de l'ordre de quelques dizaines de kHz et avec des fenêtres d'observation d'une longueur maximale de N = 64 points, c'est-à-dire d'une durée maximale de 1 à 2 ms (au lieu de 5 à 10 ms dans les appareils commercialisés actuellement). Comme on n'a pas besoin d'une analyse spectrale du signal Doppler à des intervalles de temps aussi brefs, on peut organiser les données échantillonnées en blocs de N = 64 points au maximum et ne travailler que sur un bloc sur dix. Ceci permet d'obtenir une densité spectrale de puissance toutes les 10 à 20 ms au plus, c'est-à-dire à raison d'une cinquantaine par période cardiaque.
L'initialisation peut être faite par préfenêtrage ou par covariance, comme indiqué plus haut. Le préfenêtrage, actuellement utilisé, diminue le nombre de multiplications à effectuer, mais nuit à la vraisemblance du modèle, surtout si p est proche de N. La méthode de covariance réduit l'ordre maximal du modèle à p = N/2 si on commence les récurrences à l'instant p+1 (en cas de fenêtres isolées) et elle est un peu plus complexe.
Le choix de l'hyperparamètre μ est essentiel car il règle le degré de douceur de la solution et doit maximiser la vraisemblance du modèle du signal par rapport aux données : les valeurs des paramètres AR calculés ne dépendent en effet pas des variances σ2 et t mais uniquement de leur rapport μ. L'expérience montre que la courbe de variation de la vraisemblance en fonction de μ varie très peu au voisinage du maximum et que les résultats de l'analyse sont peu affectés par des variations d'un ordre de grandeur de ce paramètre. L'expérience montre aussi que, pour un signal Doppler donné, les variations de la valeur optimale de μ d'un bloc de données à un autre sont aussi d'un ordre de grandeur. La complexité de la méthode peut donc être réduite en substituant, à l'étape de maxlmisation de la vraisemblance du modèle qui doit être effectuée sur chaque bloc de données, une étape préalable de classification des signaux en fonction des valeurs constantes à attribuer à μ. Pour que ce soit possible, il est indispensable de normaliser μ par rapport à la puissance du signal Doppler. En effet, le coefficient intervient dans le critère minimisé :
J2(a) - J1(a) + μ|a-a0|2 (7) où J1(a) a la dimension de l'énergie du signal et les paramètres AR sont sans dimension. Le coefficient μ a donc une dimension et varie comme le carré de l'échelle dans laquelle est mesuré le signal Doppler.
Pour rendre μ indépendant de la puissance du signal, il faut donc choisir : μ = μ0σ2(N/p) où μ0 est le nouveau paramètre de réglage, cette fois sans dimension, et σ2 est la variance du signal qui peut être fournie par le module d'acquisition du dispositif utilisé (24). Cette normalisation permet de travailler avec μ0 constant tout au long d'un cycle cardiaque, indépendamment des fluctuations de puissance du signal Doppler.
L'invention est susceptible de nombreuses variantes de réalisation. En particulier, on peut réduire la sensibilité du procédé à la précision limitée des calculs en propageant les racines carrées des matrices de covariance plutôt que les matrices elles-mêmes : une telle modification s'effectue de façon simple, comme indiqué plus haut. Constitution du dispositif
Le dispositif utilisé pour la mise en oeuvre du procédé peut avoir la constitution générale montrée en figure 4, qui recherche un compromis entre la qualité des résultats et le volume des calculs, est adaptable sur des appareils existants et utilise les possibilités intrinsèques de processeurs de traitement numérique du signal actuellement disponibles.
Le dispositif comporte deux processeurs couplés pour effectuer les calculs du spectre Doppler. Le premier, 20, calcule les paramètres d'un modèle du signal Doppler reçu, par l'algorithme B-CAR. Le second, 22, fournit, à partir de ce modèle, le spectre. Ce second processeur 22 commande une carte d'acquisition 24 qui acquiert les deux voies Doppler simultanément, à une fréquence d'échantillonnage qui est par exemple de 40 kHz pour l'étude de la circulation sanguine, choisie pour optimiser le compromis entre la durée des tranches analysées et la précision fréquentielle de l'analyse spectrale. Pour éviter une éventuelle dissymétrie dans l'acquisition des signaux Doppler, un seul convertisseur est utilisé. La carte peut être prévue pour détecter une éventuelle saturation des signaux Doppler.
Le processeur 20 est programmé pour mettre à l'échelle l'hyperparamètre μ de la méthode à partir de l'estimation de la puissance instantanée du signal Doppler fournie par le processeur 22, et calculer les paramètres AR du signal Doppler sur chaque bloc de données transmis par le même module.
Le processeur 22 calcule le spectre du signal sur un bloc, à partir des paramètres AR fournis par le processeur 20 sous le contrôle du module de visualisation 26 qui fixe le nombre de points sur l'axe des fréquences et le nombre de niveaux de gris ; il peut être prévu pour calculer une éventuelle composante continue dans la chaîne d'acquisition. Il doit de plus évaluer la puissance instantanée du signal, utilisée par le processeur 20 pour optimiser les résultats fournis. L'architecture de l'appareil peut être celle montrée en Figure 5. Le processeur 22 exécute l'algorithme de calcul de la densité spectrale de puissance et gère les interfaces 28-30-32 d'entrée et de sortie de données. Le processeur 20 exécute l'algorithme B-CAR.
Un mini-ordinateur hôte constitue une machine maître qui charge le programme dans les processeurs et qui commande l'exécution des programmes. Une implantation en ROM du programme permettrait d'éviter la nécessité d'un mini-ordinateur hôte. Les interfaces d'acquisition 30 et de visualisation 32 transfèrent les données, d'une part entre le dispositif d'acquisition et le processeur 22, d'autre part entre le dispositif de visualisation 26 et ce même processeur 22. Chaque processeur est muni d'une mémoire propre de travail 36.
Il est nécessaire d'effectuer les calculs en format 32 bits avec virgule flottante. Des circuits disponibles sur le marché permettent de répondre à de tels besoins. Il s'agit d'une part des processeurs de traitement du signal composés d'une unité de calcul 32 bits en virgule flottante, d'une unité de commande, de mémoires et de dispositifs d'interface intégrés sur un seul circuit (par exemple DSP32 de ATT et μPD77230 de NEC), d'autre part, des unités arithmétiques 32 bits. On a notamment implémenté le procédé à l'aide de processeurs DSP32 dont les bus internes sont accessibles de l'extérieur du boîtier (ce qui permet d' augmenter la taille mémoire de 4 Koctets à 56 Koctets), et qui ont une interface parallèle permettant d'effectuer des accès-mémoire en DMA quelle que soit l'activité du processeur et une interface série à haut débit (8 Mbits/sec).

Claims

REVENDICATIONS
1.Procédé d'analyse spectrale en temps réel de signaux instationnalres à valeurs complexes, représentant un phénomène physique, suivant lequel on modélise le phénomène représenté par le signal, par un processus auto-régressif, caractérisé en ce que :
- on modélise des tranches successives du signal échantillonné, par un processus auto-régressif d'ordre élevé, proche ou égal du maximum autorisé par la durée des tranches et la fréquence d'échantillonnage, en faisant une estimation des paramètres caractérisant le modèle par filtrage de Kalman rapide utilisant un algorithme qui est simplifié dans toute la limite autorisée par la stationnarité locale du modèle dans la tranche et en ce qu'on détermine la densité spectrale de puissance g (f) à partir des paramètres estimés ; et en ce que :
- on régularise les paramètres a du modèle auto-régressif :
P
z(n) - Σ aiy(n-i) = atyp(n)
i=1 avec at = [ap,ap-1,...., a1] yp(n)t = [y(n-p),..., y(n-1)]
où at est le vecteur des coefficients du modèle dont l'ordre est p, où b(n) représente le processus générateur du signal, où z(n) désigne les échantillons de la fonction représentant le phénomène, et où y(n) désigne les échantillons observés, par mise en oeuvre d'une contrainte consistant à minimiser un critère de la forme : J2( a) = J1( a ) + Ω( a) où J1(a) est la somme des moindres carrés [y(n)-z(n)]2, et où Ω(a) est une fonctionnelle régularisante tenant compte de la connaissance a priori de la nature du phénomène.
2. Procédé selon la revendication 1, caractérisé en ce que Ω(a) est de la forme : Ω(a) - μ |a - a0|2 où a0 est une valeur prédéterminée, choisie a priori nulle ou constituée par le résultat obtenu pour la tranche précédente.
3. Procédé selon la revendication 2, caractérisé en ce qu 'on prend μ = μo σ2 où μo est une valeur prédéterminée et σ2 la puissance estimée du signal.
4. Procédé selon la revendication 2 ou 3, caractérisé en ce que la densité spectrale de puissance est calculée par la relation : g(f) = σ2/ [ exp (-2iπfk)| 2].
5. Procédé selon la revendication 2, 3 ou 4, caractérisé en ce qu'on traite simultanément les composantes en phase et en quadrature d'un signal à bande étroite démodulé.
6. Procédé selon l'une quelconque des revendications 2 à 5, caractérisé en ce que, les signaux étant des signaux Doppler démodulés représentatifs de l'écoulement sanguin présentant un spectre de fréquence de quelques kHz, on utilise des tranches d'environ 64 échantillons prélevés à une fréquence d'environ 40 kHz et on traite simultanément les composantes en phase et en quadrature du signal démodulé.
EP90905549A 1989-03-22 1990-03-22 Procede et dispositif d'analyse spectrale en temps reel de signaux instationnaires complexes Ceased EP0464115A1 (fr)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
FR8903761A FR2644915A1 (fr) 1989-03-22 1989-03-22 Procede et dispositif d'analyse spectrale en temps reel de signaux instationnaires complexes
FR8903761 1989-03-22

Publications (1)

Publication Number Publication Date
EP0464115A1 true EP0464115A1 (fr) 1992-01-08

Family

ID=9379955

Family Applications (1)

Application Number Title Priority Date Filing Date
EP90905549A Ceased EP0464115A1 (fr) 1989-03-22 1990-03-22 Procede et dispositif d'analyse spectrale en temps reel de signaux instationnaires complexes

Country Status (4)

Country Link
US (1) US5229716A (fr)
EP (1) EP0464115A1 (fr)
FR (1) FR2644915A1 (fr)
WO (1) WO1990011494A1 (fr)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5517122A (en) * 1993-11-04 1996-05-14 The Regents Of University Of California T2 restoration and noise suppression of hybrid MR images using Wiener and linear prediction techniques
CA2331116A1 (fr) * 2001-01-15 2002-07-15 Chenomx, Inc. Identification et quantification de composes aqueux--technique et procede utilisant un systeme automatise de mesure a resonance magnetique nucleaire
US6751564B2 (en) 2002-05-28 2004-06-15 David I. Dunthorn Waveform analysis
WO2004032745A2 (fr) * 2002-10-09 2004-04-22 Ruprecht-Karls-Universität Heidelberg Calcul et representation a valeurs vectorielles de la vitesse dans un espace en trois dimensions
US20040239415A1 (en) * 2003-05-27 2004-12-02 Bishop Christopher Brent Methods of predicting power spectral density of a modulated signal and of a multi-h continuous phase modulated signal
ITFI20030254A1 (it) * 2003-10-08 2005-04-09 Actis Active Sensors S R L Metodo e dispositivo perfezionati per l'analisi spettrale
US7792314B2 (en) * 2005-04-20 2010-09-07 Mitsubishi Electric Research Laboratories, Inc. System and method for acquiring acoustic signals using doppler techniques
US8315857B2 (en) * 2005-05-27 2012-11-20 Audience, Inc. Systems and methods for audio signal analysis and modification
FR2890450B1 (fr) * 2005-09-06 2007-11-09 Thales Sa Procede de determination par analyse doppler a haute resolution du champ de vitesse d'une masse d'air
JP2013205830A (ja) * 2012-03-29 2013-10-07 Sony Corp トーン成分検出方法、トーン成分検出装置およびプログラム
US10247795B2 (en) * 2014-12-30 2019-04-02 General Electric Company Method and apparatus for non-invasive assessment of ripple cancellation filter
CN107063370B (zh) * 2017-04-18 2019-05-07 中国工程物理研究院材料研究所 一种热分布式气体质量流量计的时间响应补偿方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4391149A (en) * 1981-05-15 1983-07-05 Fischer & Porter Company Doppler-type ultrasonic flowmeter
US4585992A (en) * 1984-02-03 1986-04-29 Philips Medical Systems, Inc. NMR imaging methods
JPS6342209A (ja) * 1986-08-07 1988-02-23 Mitsubishi Electric Corp デジタル制御形自動利得制御装置
JPS6342229A (ja) * 1986-08-08 1988-02-23 Canon Inc 信号処理方法
US4789832A (en) * 1987-01-14 1988-12-06 Jeol Ltd. Three-dimensional NMR spectroscopy
US5019978A (en) * 1988-09-01 1991-05-28 Schlumberger Technology Corporation Depth determination system utilizing parameter estimation for a downhole well logging apparatus

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
See references of WO9011494A1 *

Also Published As

Publication number Publication date
WO1990011494A1 (fr) 1990-10-04
US5229716A (en) 1993-07-20
FR2644915A1 (fr) 1990-09-28

Similar Documents

Publication Publication Date Title
EP0749083B1 (fr) Procédé et dispositif de détermination du spectre de fréquence d'un signal
EP2258119B1 (fr) Procede et dispositif pour la determination de fonctions de transfert de type hrtf
WO1990011494A1 (fr) Procede et dispositif d'analyse spectrale en temps reel de signaux instationnaires complexes
EP0605281A1 (fr) Procéde de débruitage vectoriel de la parole et dispositif de mise en oeuvre
EP1615049B1 (fr) Traitement cohérent rapide pour codes à spectre de raies périodiques
EP0511095B1 (fr) Procédé et dispositif de codage-décodage d'un signal numérique
EP1410240B1 (fr) Procede et circuit d'analyse frequentielle en temps reel d'un signal non stationnaire
FR2890450A1 (fr) Procede de determination par analyse doppler a haute resolution du champ de vitesse d'une masse d'air
FR3020157A1 (fr) Procede de detection numerique
EP0166836B1 (fr) Procédé de caractérisation de la structure d'un milieu, et dispositif mettant en oeuvre ce procédé
EP3671250A1 (fr) Interféromètre numérique à sous-échantillonnage
EP1172980B1 (fr) Dispositif pour la classification de signaux complexes contenant une modulation numérique linéaire
CN112652290B (zh) 产生混响音频信号的方法及音频处理模型的训练方法
EP0673113A1 (fr) Système de caractérisation de sources de signaux
FR2549957A1 (fr) Procede et appareil pour synthetiser un signal continu estime a partir de segments d'un signal gaussien resultant de la mesure doppler par ultra-sons de l'ecoulement d'un fluide
CN113688655B (zh) 干扰信号的识别方法、装置、计算机设备和存储介质
EP4248231A1 (fr) Localisation perfectionnée d'une source acoustique
FR2854247A1 (fr) Procede de traitement de cubes sismiques correspondant pour une meme zone au sol, a differentes valeurs de deports source/recepteur et/ou d'angles d'incidence
EP3971600A1 (fr) Procédé de séparation de signaux électromagnétiques et d'estimation de leurs directions d'arrivée dans un goniomètre interférométrique large bande à réception numérique sous-échantillonnée
WO2009092881A1 (fr) Procede et dispositif pour ameliorer la resolution d'une image ultrasonore
EP0333566B1 (fr) Procédé et dispositifs de translation le long de l'axe des fréquences du module de la fonction de transfert d'un filtre
Lysenko et al. Deep learning approach to classification of acoustic signals using information features
FR2737578A1 (fr) Dispositif radar doppler a impulsions avec determination complete du vecteur vitesse de la cible
WO1990011528A1 (fr) Procede et dispositif d'extraction d'estimateurs des caracteres d'un ecoulement instationnaire a partir d'un signal doppler a valeurs complexes
EP0749014A1 (fr) Procédé de reconstruction d'un champ dense de caractéristiques associées à un phénomène physique

Legal Events

Date Code Title Description
PUAI Public reference made under article 153(3) epc to a published international application that has entered the european phase

Free format text: ORIGINAL CODE: 0009012

17P Request for examination filed

Effective date: 19911011

AK Designated contracting states

Kind code of ref document: A1

Designated state(s): DE DK ES FR GB IT NL

17Q First examination report despatched

Effective date: 19930305

STAA Information on the status of an ep patent application or granted ep patent

Free format text: STATUS: THE APPLICATION HAS BEEN REFUSED

18R Application refused

Effective date: 19930827