EP0422158A1 - Process and device for extracting certain parameters from the features of an unsteady flow based on a complex-value doppler signal - Google Patents

Abstract

The method and the device can in particular be used to extract the estimators of the Doppler signal in ultrasound medical imaging. The device uses as input parameter, to calculate the average frequency, the correlation coefficients Im (Cxx) corresponding to the imaginary part of the signal. These coefficients are applied in sequence to an input of a multiplier (10). The other entry receives weights of the form 0, 1/0, 1 / pi, ..., (- 1) p + 1 / pip, where p is the number of weights, stored in a register (11). The partial products are accumulated and provide a sum representative of the average frequency.

Description

The invention relates to the extraction of character estimators of an unsteady flow from a Doppler signal with complex values. The invention relates to the extraction of character estimators of an unsteady flow, from a signal. Sampled complex value Doppler. The term "estimator" designates a parameter providing an indication representative of a characteristic of the flow at a determined time, such as in particular the average Doppler frequency (representing the average speed of the flow) and the average bandwidth of the spectrum. Doppler (representing the dispersion of velocities in the flow).

 The invention finds a particularly important application in medical ultrasound imaging with Doppler effect: the average frequency then gives an indication of the speed in a vessel or a cavity. The bandwidth makes it possible to determine the flow regime and, in particular, to detect the jet regimes.

There are calculation methods making it possible to determine the average frequency and the average bandwidth of a Doppler signal, for a determined block of data on which the signal is almost stationary, from the power spectral density or dsp, g _Y (f) for the data block.

The average frequency is given by

(1) The average bandwidth Δf is given by (2)

 'J

where is the average frequency and

the average of the squares of given by

 /

(2bis)

 Formulas (2bis) and (1) assume that the d.s.p. is previously normalized, i.e. we have

We know a method of calculating the average frequency by an autocorrelation technique (Kasai et al, "Real time two-dimensional blood flow imaging using an auto-correlation technique", IEEE TRANSACTIONS on sonics and ultrasonics, Vol. 32 , pp. 458-464, May 1985). It has the disadvantage of starting from the assumption that the modulus of the correlation function of the signal considered has a zero derivative at the zero point, whereas this assumption is far from being always verified.

Furthermore, patent application FR n ° 89 03761 for "Method and device for real-time spectral analysis of complex unsteady signals" makes it possible to represent the signals concerned by the present invention using an autoregressive parametric model, whose coefficients remain stationary over short time intervals. We could therefore calculate the estimators mentioned above from the coefficients of the autoregressive model. But this approach results in a large volume of computation and difficulties of computation of integral.

 The invention aims to provide a method for extracting estimators making it possible to arrive at high precision by calculations on a small number of samples, therefore with a relatively reduced calculation volume for a significant estimator.

To this end, the invention uses the covariance function of the sampled signal as a calculation aid. The latter being locally stationary, its covariance is stationary in each block of samples and identifies with an autocorrelation function whose coefficients C _Y (n), for whole delays, are the Fourier coefficients of the dsp g _y (f) which is itself a periodic function of the frequency. So we have :

in which :

We see that we replace the integrals (1) and (2bis) by weighted sums of the coefficients of the signal autocorrelation function. In formula (4), the weight assigned to successive values C _Y (n) of the correlation function decreases very quickly with n, so that a small number of terms may suffice, when the coefficients C _Y have the same order of greatness.

We thus arrive at a value of the average frequency given by:

We see that the average frequency can be calculated in a simple way from the coefficients of the correlation function.

 These coefficients can be calculated:

- either directly from the samples,

- either from the coefficients of an autoregressive parametric model, by a fast algorithm, for example of the "inverse Levinson" type which will be recalled below.

 In the latter case, it is possible to use, without the drawbacks mentioned above, the results provided in the patent application of the same day already mentioned.

 The invention proposes a method for extracting estimators according to which: the Doppler spectrum of a flow is measured and the Doppler signal is sampled, characterized in that: the coefficients of the autocorrelation function are calculated for each block d 'samples; the sum of the correlation coefficients corresponding to the different available delays is calculated, each weighted by the integral of the product of the frequency by a negative exponential as a function of the delay to obtain the average frequency.

We can consider the average frequency as the 1st order moment of the spectrum of the complex Doppler signal, the bandwidth of the spectrum as 2nd order moment of the spectrum and other characteristics (asymmetry of the spectrum or flattening) as represented by higher order moments.

It is possible to significantly simplify the extraction calculations by analyzing and adapting formula (5). This formula shows that the calculation of the average frequency boils down to calculating a sum of correlation coefficients weighted by integrals of the form:

(6)

Now, for n = 0, we have an integral of the form f.df which is zero and we can write, in general, the integrals for n 0 in the form:

(7) We deduce from this, taking into account the properties of Hermitian symmetry of the autocorrelation coefficients of a signal with complex values, formula (8):

(8) We see that it is enough, to calculate the average frequency ee have the imaginary part Im (Cxx) autocorrelation coefficients of the sampled input function.

 We can also notice that formula (8) shows a real average frequency, which corresponds to physical reality and that, in the case of a signal with real values, whose correlation function is real and even, the frequency mean is zero.

 The invention consequently also proposes a method according to which the sum of the imaginary parts of the autocorrelation coefficients of the signal is each assigned a weighting coefficient.

 If the number of available coefficients were infinite, the value of the weighting coefficient would only be a function of n for the best possible approximation.

However, in practice, the number p of coefficients is limited to a low value which leads to deviations for the frequencies close to the sampling frequency. Since it is often sufficient to obtain a satisfactory representation in a limited frequency range, where the dsp is maximum, the factors 1 / n can be replaced by factors a _n chosen to improve the representation in the zone judged essential, at the cost of degradation in other areas. The selection of a _{n will be} carried out by a process of optimization of a polynomial approximation of the part considered essential.

 The average frequency is thus obtained directly.

The invention also provides a device for implementing the above-defined method, which can be entirely wired and therefore capable of operating in real time at high speed, comprising only elements of a simple nature (multipliers and accumulators). The invention will be better understood on reading the following description of particular modes of implementation, given by way of nonlimiting examples. The description refers to the accompanying drawings, in which:

 - Figure 1 is a block diagram of a circuit for calculating the average frequency and the average bandwidth of an unsteady Doppler signal with complex value, representing the Doppler spectrum of a flow from the real parts and imaginary of its correlation function;

 - Figure 2 is a block diagram showing a circuit for calculating the imaginary part of the coefficients of the signal autocorrelation function, coefficients acting as inputs for the calculation of the average value in the circuit of Figure 1;

 - Figure 3, similar to Figure 2, shows a circuit for calculating the real part of the coefficients of the autocorrelation function, involved in the calculation of the average bandwidth.

We will successively describe the calculation of the average frequency and that of the average bandwidth Δf using the device according to the invention shown in Figure 1. It should however be noted that, often, it is not necessary for the user of the device to have both the average frequency and the bandwidth and that, in this case, the device can be simplified by removing unnecessary components.

Calculation of the average frequency

In the embodiment of Figure 1, the average frequency is calculated by bringing into play all the terms of formula (8). But it will often be possible to truncate the sum by omitting the terms corresponding to a value of n greater than a determined threshold. Indeed, the weight of successive correlation coefficients varies in 1 / n. These coefficients themselves cannot have a module greater than unity and are often of the same order of magnitude. Consequently, their contribution often decreases rapidly with n.

In Figure 1, the notation Im (Cxx) designates the coefficients Im which are involved in formula (8). To allow synchronization of the entire circuit (from a central clock not shown), the successive coefficients pass through a cascade constituted by a memory S and a network of flip-flops or "latch" L. The coefficients corresponding to the delays successive are applied one after the other to a multiplier 10 whose second input receives the multiplier coefficients 0, 1 / π, ... (-1) ^{n + 1} / πn stored in a register 11 or an addressable memory for form the successive terms of the sum of formula (8).

 The successive products are accumulated for all the values of n = 1 to n = + p in an adder constituted for example by a shift register 12 looped by a network of rockers 14 forming "latch".

 The average frequency appears in a memory

16 at the exit of the register. According to formula (8), this value is normalized as a function of the signal strength of which the auto-correlation moment for a zero delay constitutes a representation. In Figure 1, this zero delay autocorrelation moment is also used for the calculation of the average bandwidth. The autocorrelation coefficients of the signal are taken directly from the input of the real parts Re. Normalization can be carried out using a divider or, better, as indicated in FIG. 1, using a table 18 constituted by a logical network programmable.

 This fundamental process can be refined to further improve the quality of the estimate of the average frequency, by acting on the estimator and / or the iterations of the estimator.

 I - The estimator can be modified to correct various effects and improve the compromise between the robustness with respect to the approximations introduced by the windowing of the signal and the frequency dynamics.

1. Correcting the windowing effect on the correlation means amounts to correcting the bias introduced into the calculation of the correlation coefficients by the windowing of the signal. The calculation performed is not that of Cy (n) but that of ((pn -) / p) Cy (n), where: p is the number of points in the window,

n the delay in the correlation coefficient. The effect of Cy (1) being by far the most important in the estimation of f, we propose as estimator: f ₁ = (P / (P-1)) f

f _χ = (P / (P-1)) ∑ [(1 / πn) (-1) ^{n + 1} Im <Cy (n)>]

The summation is performed for n varying from 1 to p.

2. The correction of the windowing effect on the weighted sum of the correlation coefficients is carried out by using more correlation moments in the calculation of the mean frequency.

For a better estimate, we replace the sum: f = Σ [(1 / πn) (-1) ^{n + 1} Im <Cy (n)>)] with the summation performed for n varying from 1 to p by the sum: f ₂ - ∑ [(1 / π) (-1) ^{3m-1 + 1} aα (3m-1) im <C'y (3m-1) >].

The summation is carried out for m varying from 1 to p and 1 varying from 0 to 2.

 The values C'y (m-1) are obtained by non-linear extrapolation from the values Cy (n).

 In what follows, the correlation function

Cy (n) will be considered as a complex function of the form:

R (n) exp (iθ (n)) of module R (n) and of argument θ.

R (n) = [(Re <Cy (n)>) ² + (Im <Cy (n)>) ² ] ^1/2

R (n-1) = [(Re <Cy (n-1)>) ² + (Im <Cy (n-1)>) ² ] ^1/2

 R (3m) = R (n)

 R (3m-1) = [2R (n) + R (n-1)] / 3

R (3m-2) = [R (n) + 2R (n-1)] / 3 θ (n) = tg ^-1 [(Im <Cy (n)>) / (Re <Cy (n)>) ]

θ (n-1) = tg ^-1 [(Im <Cy (n-1)>) / (Re <Cy (n-1)>)]

θ (3m) = θ (n)

θ (3m-1) = [2θ (n) + θ (n-1)] / 3

θ (3M-2) = [θ (n) + 2θ (n-1)] / 3

C'y (3m) = R (3m) [cosθ (3m) + i.sinθ (3m)]

 C'y (3m-1) = R (3m-1) [cosθ (3m-1) + i.sinθ (3m-1)]

C'y (3m-2) = R (3m-2) [cosθ (3m-2) + i.sinθ (3m-2)].

The coefficients aα (m-1) are calculated as follows, for real α between [0 and 1] and set by the user. aα (3m) = 3 / 3m = 1 / m

aα (3m-1) = 3α / (3m-1)

aα (3m-2) = 3α / (3m-2). Note that if α = 0: f ₂ = Σ (1 / π) [(-1) ^{3m + 1} (1 / 3m) Im <C'y (3m)>]

f ₂ = f because n - 3m. 3. Finally, the correction of the windowing effect on the correlation moments and on the weighted sum of the correlation coefficients is carried out by cumulating the two preceding improvements 1 and 2: f ₃ = [(1 + 2α) p / ((1 + 2α) p-1)] f ₂

f ₃ = [(1 + 2α) p / ((1 + 2α) p-1)] x Σ [(1 / π) (- 1) ^{3m-1 + 1} aα (3m-1) Im <C'y (3m-1)>]. It should be noted that, if α = 0: f ₃ = f ₁

If a = 1 f ₃ _ (1 / 2π) tg ^-1 [Im <Cy (1)> / Re <Cy (1)>].

II - The iterations of the estimator can be modified to overcome the Robustness / Frequency dynamics compromise of the previous estimators.

1. At the level of local estimation, an iterative shift of the integration domain is used in the calculation of the average frequency so as to overcome the Gibbs phenomenon linked to the truncation of the sum over the correlation moments. We have p samples of the signal isolated as above.

We use one of the four previous estimators f, f ₁ , f ₂ or f ₃ called f _u in the following.

 We loop the estimator on itself according to the following sequence:

Initialization:

k = 0

f _u 0 = 0

Index increment:

k = k + 1 Frequency shift of the Doppler signal:

y ^{k + 1} (t) = y ^k (t) exp (-2πf _u ^k t)

Calculation of the Doppler signal correlation moments: Cy (n) or C'y (m-1)

Calculation of the average frequency:

f _u ^{k + 1} = f ^{k + 1} , f ₁ ^{k + 1} , f ₂ ^{k + 1} or f ₃ ^{k + 1}

End of iteration test:

Sign [f _u ^{k + 1} ] 4 Sign [f _u ^k ] or:

k> k ₀

k ₀ defined in advance by the user, or:

the two tests combined.

 If the test is validated, then:

f = f _u ^{k + 1}

Otherwise, return to incrementing the index. 2. At the dynamic estimation level, the average frequency already estimated at the previous instant, at the previous pixel of the same line or finally at the previous line of the same image is used to calculate the average frequency.

f _u ^k represents one of the estimators f, f ₁ , f ₂ or fg at time k or at pixel k or at line k.

 (a) To improve the frequency "smoothness", the estimated average frequency is not to vary too much between k and k + 1.

 We use for this the following synoptic:

Initialization:

k = 0

f _u 0 = 0

α ₀ , α ₁ , a ₂

α = α + α ₀ α ₀ , α ₁ and α ₂ are set by the user.

α ₀ is the increment used to modify α.

α ₁ prohibits rapid variations of the average frequency around zero.

α ₂ regulates the softness constraint. Increment of the window, pixel or line index: k = k + 1

Calculation of the value of o: α = α-α ₀ . Calculation of the Doppler signal correlation coefficients: C'y (m-1)

Calculation of the average frequency: f _u ^{k + 1} = f ₂ ^{k + 1} or f ₃ ^{k + 1} End of iteration test: abs [(f _u ^{k + 1} -f _u ^k ) / (f _u ^k + 1l )] <α ₂ or α <0 If the test is validated: f = f _u ^{k + 1}

Otherwise, return to incrementing the index. It should be noted that, if α ₀ = 1, the sequence of operations is reduced as indicated above in 1.3.

Estimated average frequency by the estimator: f ₄ = (1 / 2π) tg ^-1 [Im <Cy (1)> / Re <Cy (1)>]

Test: abs [(f _u ^{k + 1} -f _u ^k ) / (f _u ^k + α1)] <α ₂ . If the test is validated, we keep f.

Otherwise, the average frequency is calculated by the estimator f ₁ .

(b) To improve frequency monitoring, the average frequency at time k + 1 is calculated as the sum of the average frequency at time k and a change in average frequency at time k + 1 according to l '' following algorithm: Initialization:

k = 0

f _u ⁰ = 0

Increment of window, pixel or line index:

k = k + 1

Frequency shift of the Doppler signal:

y ^{k + 1} (t) = y ^k (t) exp (-2πf _u ^k t)

Calculation of the Doppler signal correlation moments: Cy (n) or C 'y (m-1)

Calculation of the incremental average frequency:

f _u ^{k + 1 =} f ^{k + 1} f ₁ ^{k + 1} f ₂ ^{k + 1} , or f ₃ ^{k + 1}

Calculation of the average frequency

f = f _u ^{k + 1} + f _u ^k

Calculation of average bandwidth

The average bandwidth used to calculate the circuit illustrated in Figure 1 is defined by formula (2) for a standardized dsp. The average frequency being known, it remains to calculate (or the variance σ) which constitutes the second order moment of the dsp:

(9) who can write

(10)

The calculation is very similar to that of : it consists in making the sum of autocorrelation coefficients C _Y (n) weighted by the integrals:

(11) The integral is equal to 1/12 for n = 0; for n 0, the integral of order n is equal to:

(12)

Still using the Hermitian symmetry properties of the auto-correlation functions, we arrive at expression (13) similar to formula (8):

(13) The terms 1 / n ² can be replaced by terms b _n ² to ensure an optimal approximation in a particularly interesting area (at maximum dsp).

We see that the mean of the square of the frequency reduces to the constant term 1/12 if the signal is white noise, which is normal since the dsp remains equal to unity over the interval [-1/2, + 1/2]; the average frequency is then zero, Δf takes its maximum value which is equal to 1 /

The part of the circuit of Figure 1 providing the mean of the square of the frequency has a constitution similar to that of the part of calculation of The only difference is to add 1/12 to each input value from the looped register 12a, before normalization.

The calculation of Δf from the available values of and ² can be carried out in a simple manner in a circuit 20.

As indicated above, the average value can be calculated from only the imaginary parts of the autocorrelation coefficients which can be calculated directly from the samples of the signal, which has a real part Re and an imaginary part Im.

 The circuit shown in Figure 2 has two channels with the same constitution. The samples Re and Im of the Doppler signal are applied to two respective inputs 22 and 24. To allow the circuit to process blocks of successive and adjacent samples, each channel comprises, in FIG. 2, two sets of shift registers on each track, operating alternately. The samples Re of the real part are alternately directed, by the switch 30, to two registers 26 and 26a and two registers 28 and 28a. Similarly, the samples Im are alternately directed, by a switch 32 controlled in synchronism with the switch 30, to registers 34 and 34a and registers 36 and 36a.

 Multipliers 38 and 40 receive the signals from the registers through switches controlled in synchronism with switches 30 and 32 and carry out the partial sums of the samples of the same block, designated by v1 to vn in Figure 2, are first in phase, then with successive shifts.

 The diagram shown in Figure 2 does not include means stopping the calculation when the terms become negligible; the calculation is therefore continued with n successive shifts, in a number equal to that of the available samples.

 The calculation of the successive terms Im (Cxx) of the correlation function is carried out at each point using a subtractor 42; the terms are accumulated in a register 44 looped by a network of flip-flops 46 until the totalization has been carried out on the n points available.

 We thus find successively, on the output Im

(Cxx) the imaginary part of the first moment (after the first calculation sequence) at the end of the first sequence, then the imaginary parts of the other moments, each time corresponding to a determined offset, during the following sequences. These moments are successively applied, in the diagram of Figure 1, to the multiplier 10.

 If it is necessary to calculate the real autocorrelation coefficients, we can use a circuit of the kind shown in Figure 3, which will not be described given its similarity to that of Figure 2.

 The circuit shown in Figure 1 can be completed to go beyond the second order moment, in particular in the case of a highly disturbed signal or with a very asymmetrical spectrum.

 An alternative embodiment of the invention uses, as input data, not the samples of the real and imaginary parts of the complex signal, but the parameters of an autoregressive model obtained for example by implementing the method according to the request âe Patent "Method and device for real time spectral analysis of complex unsteady signals" FR n ° 89 03761.

For this, it implements a Levinson algorithm which makes it possible to calculate the coefficients constituting the correlation matrix from the matrix of the parameters of the autoregressive model of the Doppler signal. A description of Levinson's algorithm can be found in G. Demoment, Algorithmes Rapides, Cours de l'Ecole Supérieure d'Electricité, n ° 3152, 1987. The inverse Levinson algorithm is a fast algorithm requiring order of p ² elementary arithmetic operations to pass from a prediction vector a _p to the corresponding correlation matrix. The matrix obtained is a Toeplitz matrix satisfying the optimality conditions and guaranteeing high precision, given the absence of approximation. The application of the method described above to medical Doppler imaging has shown satisfactory results on short blocks, of only eight samples, representative of current working conditions in color Doppler imaging. In particular, the method according to the invention revealed a lower variance than that of the conventional methods for determining the average frequency by correlation methods, on the estimation of the average frequency in diastole.

Claims

1. Method for extracting character estimators, and in particular the speed of an unsteady flow, according to which the Doppler spectrum of a flow is measured and the Doppler signal is sampled, characterized in that: the coefficients are calculated the autocorrelation function for each block of samples; the sum of the correlation coefficients corresponding to the different available delays is calculated, each weighted by the integral of the product of the frequency by a negative exponential as a function of the delay to obtain the average frequency.

 2. Method according to claim 1, characterized in that the sum of the imaginary parts of the signal correlation coefficients is calculated, for each correlation coefficient, a weighting coefficient whose value is a function of the order of autocorrelation coefficient.

 3. Method according to claim 1 or 2, characterized in that the bandwidth is determined by calculating the average of the squares of the frequency and subtraction.

4. Device for extracting estimators of the characters of an unsteady flow, from a sampled Doppler signal with complex values, characterized in that it comprises: means for applying in sequence the imaginary parts of the correlation coefficients of the signal during each time said interval, at an input of a multiplier (10); means (11) for applying to the other input of the multiplier (10), in synchronism, weighting coefficients of the form 0, 1 / π, ..., (-1) ^{P + 1} / πp, where p is the number of weights; and means (12, 14) for accumulating the partial products and providing a sum representative of the frequency average.

 5. Device according to claim 4, characterized in that it further comprises means (18) for normalizing the average frequency, comprising a transcoding table (18) which receives the output of the accumulation means and a value representative of the autocorrelation coefficient of the function for a zero delay.

6. Device according to claim 4 or 5, characterized in that it further comprises means for determining the average of the squares of the frequency, comprising means for applying in sequence the real parts of the correlation coefficients of the signal during each times said time interval at an input of a multiplier; means for applying to the other input of the multiplier, in synchronism, weighting coefficients of the form 0, -1 / π ² , ..., (-1) P / p ² π2; and means for accumulating the partial products.

 7. Device according to claim 6, characterized in that it further comprises means for normalizing the average of the squares of the frequency, comprising a transcoding table (18a) receiving the autocorrelation coefficient for a zero and expected delay to add a constant value to the output of the accumulation means and divide the result by said zero delay autocorrelation coefficient.

 8. Device according to any one of claims 4 to 7, characterized in that it further comprises means for directly determining the imaginary parts of the correlation coefficients, and possibly the real parts of the correlation coefficients, from samples of the real part and the imaginary part of the signal, each time over said time interval.

9. Device according to any one of Claims 4 to 7, characterized in that it further comprises means for determining said correlation coefficient from the parameters of an autoregressive model of the signal by application of the inverse Levinson algorithm.

 10. Device according to claim 4 or 6, characterized in that the weighting coefficients are adjusted to optimize the approximation of the average frequency and / or the bandwidth in a determined and limited frequency range.