- Email: [email protected]
Adaptive estimation of the mean frequency of a doppler signal from short data windows
Ultrasoundin Med. &Biol. Vol. 17, No. 9, pp. 901-919, 1991
0301-5629/91 $3.00 + .00 © 1991 Pergamon Press plc
Printed in the U.S.A.
OOriginal Contribution ADAPTIVE ESTIMATION OF THE M E A N FREQUENCY OF A DOPPLER SIGNAL FROM SHORT DATA W I N D O W S A. HERMENT, t G. DEMOMENT,* J. P. GUGLIELMI, t PH. DUMEE t a n d C. PELLOTt tINSERM U.256, H6pital Broussais, Paris, France; *L2S(CNRS/ESE), UPS, Gif-sur-Yvette, France (Received 7 March 1991; in final form 8 August 1991 )
Abstract--The color Doppler estimator ( C E 1), which is calculated from the phase of the first correlation lag of the Doppler signal, is compared to the general mean frequency estimator (CE.), which is based on a weighted summation of all the available correlation lags, for long and short Doppler data sets (typically 48 and 8 Doppler samples). A new estimator of the Doppler signal mean frequency is derived from the results of this study. It optimizes the compromise between the range of analyzable frequencies and the estimation variance for the characteristics of the Doppler signal. Demonstration is provided that the behavior of this estimator shifts from that of CEI to that of CE,, according to the setting of a single parameter. An adaptive version of this estimator is implemented and applied to Doppler recordings. Applications can be contemplated for color Doppler imaging.
Key Words: Color Doppler, Mean frequency, Correlation function.
INTRODUCTION
ties of the correlation function of the Doppler signal, in particular, the double correlator (Gerzberg and Meindl 1980a) and the single correlator (Arts and Roevros 1972). This article presents a comparison between two of these correlation-based estimators of the Doppler signal mean frequency. The first estimator is the conventional color Doppler estimator (CE~), introduced by Miller and Rochwarger (1972), and also referred to as the correlation angle estimator, which utilizes the phase of the first correlation lag of the Doppler signal. The second estimator is the general mean frequency estimator (CE,), described by Angelsen and Kristoffersen (1983a), and also referred to as the mean frequency estimator, which utilizes a weighted sum of all the available correlation lags of the Doppler signal. This study will be used as a starting point to new adaptive mean frequency estimator of the Doppler signal. First, the mathematical background of these two estimators and its inference on their respective behavior for long data sets are presented. Second, a simulation study is used to quantify their performance when applied to a reduced number of Doppler data corrupted by noise. A survey of possible improvements of the general mean frequency estimator is then provided, from which an adaptive mean frequency estimator is derived. Finally, application of this estimator to actual Doppler recordings is given.
Since the early stages of ultrasound velocimetry, numerous time-domain estimators of the mean frequency of the Doppler signal have been developed to quantify the mean velocity of the flow. The first estimators were implemented on analogue systems using zero-crossing counters (Rice 1963), (Angelsen and Kristoffersen 1983b) or phase locked loop techniques (Sainz et al. 1976), (Bramanti and Marchesini 1983 ). Thereafter, other signal phase detectors such as the infinite gate detector (Nowicki and Reid 1982), the I/Q algorithm (Barber et al. 1985 ), the double-correlation centroid detector (Reid 1966 ), the single-correlation centroid detector (Arts and Roevros 1972) and, finally, the ~ f centroid detector (Gerzberg and Meindl 1980b) have been used. Recently, estimators with better, but computationally intensive, performances have been developed using autoregressive signal modelling (Kuc and Li 1985); (Loupas and McDicken 1990) or autoregressive moving average modelling of the Doppler signal (Talhami and Kitney 1988). Finally, estimators derived from the properties of the auto-correlation function of the Doppler signal have been developed and applied to color Doppler imaging. Their main advantage consists in offering the best compromise between the quality of mean frequency estimation and the shortness of available signal windows. It must be pointed out that most of the previous estimators can be derived from the proper901
902
Ultrasound in Medicineand Biology
METHODS
and, in particular, for r = 0:
Mathematical background Both estimators CE~ and CE, utilize the properties of the correlation function of a complex value signal. However, the derivation of each mean frequency estimator is based on different assumptions concerning the nature of the signal. Color Doppler estimator CEt The Doppler signal is assumed to be a complex value and continuous time signal: z(t) = x(t) + iy(t)
tE R
( 1)
where R represents the ensemble of real numbers, x(t) and y(t) are the in-phase and quadrature components (I and Q) of the analytic signal, resulting from the baseband conversion of the backscattered radiofrequency signal. This signal can either be assumed deterministic with finite energy, or random with second-order stationarity; subsequently, this difference will have no importance. In all cases, its energy spectrum (deterministic context) or its power spectral density (PSD) (stochastic context), gz(f), is a positive, real value function. The mean frequency of the Doppler signal is defined by:
f=
f
-oD
fg.(f)df
provided that its energy, in a deterministic context, or its mean instantaneous power, in a stochastic context, are normalized:
L
oo
g~(f)df = 1
(3)
~ fg~(f)df
(6)
Therefore, according to eqn (2):
l [d<(,) 1 tTJ,_o
/= E
(7)
Miller and Rochwarger (1972) and Kasai et al. (1985) utilize a modulus and phase representation of
Cz(~): (8)
G ( r ) = I G ( r ) l e '~*) Equation (7) can then be rewritten as follows:
l [ dlG(r)lei't" 1 l[[dlCz(r)l i d~(r) ] ei~(,)] - 2i~ [ -~r + ilG(r) dr J j,=0 (9) Because the magnitude of the Doppler signal IG(OI is an even function, its derivative d] Cz(r) I vanishes dr at r = 0, then:
f:27l[
1
IG(,-)I
l [ ,oU)l
d~ ],=o=2--~[
dr
],=o (lO)
Because [ cf. eqn (3) ]: IC(r)l,=o =
f
--oo
g~(f)df = 1
(11)
o()
of)
Here, E represents the expectation operator. The mean frequency can be determined from the autocorrelation function of the signal, which is the inverse Fourier transform of the PSD:
C=(r) =
= 2i~r
TJ,=o
(2)
c¢3
E{ Iz(t)l 2} =
Volume17, Number9, 1991
g,(f)e2i*Y'df
(4)
Because the phase is an odd function, ¢(r) vanishes at r = 0 so that:
¢(~) lim~o - -
(12)
1 Im{ Cz(z) } lim,-otg- Re{Cz(r) }
(13)
f =~
1
which finally leads to:
c~
Indeed, if the first time-derivative of Cz(r) does exist, we have:
dCz(r) f f o~ d ~ - 2i7r f g z ( f ) eZ"~df oo
(5)
f= ~
l
In this expression, tg -~ represents the arc tangent function. At this stage of the demonstration, one can observe that the C E 1 estimator is perfectly defined for
Adaptive estimation for Doppler signal • A. HERMENTel aL all f G [ - oo, + oo ]. There is no limitation along the range of analyzable frequencies. But, practically, the mean frequency computation must be achieved from discrete samples of the Doppler signal. One can only estimate the correlation function at discrete time instants n = 0, _+1, +_2 • • .. The mean frequency is therefore approximated
903
+oo ~+1/2 ~ C~(n) fe-Zi'f"df n=-oo ~' -1/2
f=
The mean frequency is thus expressed as the sum of the different correlation lags C~(n) weighted by the following coefficients:
by:
+1/2
f f = ~(1)
1
Im{ C~(I)}
(19)
fe-2"f"df=
-1/2
(14)
--
2---~ - 27r tg-I R e { C z ( l ) }
0 for n = 0 (-1).+1 (20) 2i~rn f o r n 4 =0
_
_
which gives: This new expression of f will induce a number of constraints on CEI, which are directly linked to signal sampling. In particular, the range of analyzable frequencies will be reduced to: [ - P R F / 2 , + P R F / 2 ] , where PRF represents the pulse repetition frequency of the Doppler system. Indeed, this first order approximation of the derivative would only be valid provided that the phase is relatively smooth or that it does not change faster than the PRF. Hence, a striking fact about CEI is that the theoretical derivation is based on a continuous signal model, and that eqn (14) is the consequence of the discrete nature of a pulsed Doppler system.
Here, the Doppler signal is assumed to be a complex value and a discrete time signal:
nEZ
(15)
where Z represents the ensemble of integer numbers. The Doppler signal is considered as a sampled version of a continuous time and band-limited signal. Whatever the case (deterministic or random signal), its PSD is a positive, periodic, real value function that can therefore be expanded in a Fourier series:
gz(f) = ~ Cz(n)e -2i'fn
(16)
nEZ
whose Fourier coefficients are defined as:
+l/2
C~(n) =
f -1/2 g=(f)e2~'/"df
nE Z
__
2i7r
X [ ~ -(-l)n+1 Cz(n) + ~ tl~-
oo
.=1
n
(-1)"+1 ] n C=(n) (21)
Taking into account the hermitian symmetry property of the correlation lags of a complex value process:
C=(n) = C*(-n)
(22)
we have:
General mean frequency estimator of CE,
z(n)=x(n)+iy(n)
1 f=
(17)
These coefficients are the correlation lags of the discrete correlation function of the signal. Because g=(f) is periodic, the mean frequency is defined differently, and eqn (2) must be rewritten as follows:
1 ~ ( - 1 ) "+t - f = 25g n
[ C z ( n ) - Cz*(n)]
(23)
and finally: f=l
~ (-l)n+______ ~ Im{Cz(n)}
"/I"n=1
n
(24)
This is an expression of the estimator, previously given by Angelsen and Kristoffersen (1983a). This expression can also be derived from the derivative of the correlation function expressed as a Taylor series expansion. But in spite of this common starting point, the different assumptions made on the nature of the Doppler signal and the way in which the derivative at the origin of C=(T) is estimated in the continuous case, give rise to very different properties for CEI and CE,.
Differences between CEI and CE. Ifa phase and modulus representation is used as in eqn (8), and if N represents the number of available correlation lags, (24) will give:
+zo Cz(n)e-2i'f"df f = ffl5~ fgz(f)df = ~'f+1/2--1/2f,,=~-oo (18)
f=l
~ (-l)"+_____~im{iCz(n)[e,~.)} (25)
71"n=l
n
904
Ultrasound in Medicine and Biology
or;
Volume 17, Number 9, 1991 Cz(n)
1 --
y=
l -
~
(-l)
71" n = l
"+'
IC~(n)l .sin ~o(n)
(26)
1
f = ~ ~o(1)
(27)
Equation (26) will reduce to eqn (27) when: N = 1, ICz(l)l =½
andsin¢(l)=~o(l)
(28)
which is very restrictive. CE~ and CE, are quite different estimators, even if the basic principle of their derivation is the same.
Influence of the correlation lag estimator At this stage the estimators are not completely defined. In practice, the correlation lags are unknown and must be estimated from the Doppler signal samples. The way these correlation lags are computed will influence the properties of the estimator. There are numerous approaches to calculate the correlation lags. However, because of computation time requirements, only two correlation function estimators can be used: 1 N-n
C,(n) ~ E Z,+mZm m=l
forn=0, I...N-
1 (29)
+PRF/2
* Zn+mZm
1 (30)
The expression z* represents the complex conjugate of z. C1 is the conventional correlation estimator that is biased but whose variance is acceptable. The C2 estimator is unbiased. However, its variance becomes high when n is large and close to N. The ratio Im{ C( 1) }/ Re {(C( 1 ) } does not depend on the type of correlation estimator; C~ or C2 can be equally used for CE~. On the other hand, selecting either Ci or C2 for CEn must still be defined. C~ will reduce the contribution oflm{ Cz(n) } with n, thus allowing more rapid convergence of the sum ofeqn (24) and a reduction of the Gibbs phenomenon associated with the truncation. However, this decreasing contribution of Im { (C~(n)} will lead to a bias of the frequency characteristic of CE,, increasing toward the high frequencies. This phenomenon is represented in Fig. 1, which displays the behavior of CE,, for a sine signal with a frequency varying from - P R F / 2 to +PRF/2 and a N = 48 sample window, when the correlation estimators Cl and C2, given in eqns (29) and (30), respectively, are applied. Because CE~ is not dependent on the exact nature of the correlation function estimator and because of the oscillations of the frequency characteristic of CEn associated with the use of the correlation estimator C2, we will subsequently adopt, in this article, the correlation estimator Ci for both the mean frequency estimators. By making this choice, no attention has been paid to the statistical properties (bias and vari-
+PRF/2
J
-PRF/2
Z
forn=0,1..-N-
n
while eqn (14) gives the estimator CEv
N-n
N-n,n= 1
A
i
i
I
+PRF/2
i
-PRF/2
B
•
•
+PRF/2
Fig. 1. Influence of the correlation estimator on CE,,, sine signal, N = 48: x = actual Doppler frequency, y = estimated Doppler frequency. (A) Cj: smooth frequency characteristic, (B) C2 oscillating frequency characteristic.
Adaptive estimation for Doppler signal • A. HERMENTel al.
ance) of the correlation lag estimator. However, the low variance of C~ will be in accordance with the statistical and noisy nature of the Doppler signal.
Theoretical behavior of the estimators An extensive analysis of CEI and CEn has been made by Kristoffersen (1988). The derivation of analytic expressions for the bias and the variance of these estimators in a hypothesis of strong filtering led to the following: CE1 yields a low bias estimate of the mean frequency, which becomes unbiased for symmetrical Doppler spectra. However, the variance of estimation increases rapidly with Doppler signal bandwidth and noise. Finally, these quantities are shown to be independent of the signal mean frequency as far as no aliasing occurs. CE~ is characterized by a bias increasing with the spectrum skewness, the noise and the Doppler signal bandwidth. This bias can reach up to 50% for very large Doppler spectra. However, the noise immunity of the estimator is far better than the one of CEi and, thus, even more when the signal bandwidth is large. Finally, the estimation variance is mean frequency independent as far no aliasing occurs. In addition, Kristoffersen and Angelsen (1985) propose an alternate approach that consists in shifting the spectrum down to the zero frequency by an adaptive premixing of the Doppler signal. This provides an unbiased mean frequency estimator with a variance that is independent of the mean frequency. However, the premixing operation can be applied in a safe way only if the strong filtering hypothesis is valid. This will not always be the case for noisy, wide-band and short data sets of Doppler signal. This is why this adaptive version of CE~ will not be studied in the following developments. However, its performances will be compared to those of other mean frequency estimators on actual Doppler signal recordings. These theoretical results (Kristoffersen 1988) are now illustrated on a synthetic Doppler signal for further comparison with short data sets.
Behavior of the estimators on long data sets Two characteristics are of major importance for a good estimation of the mean frequency of the Doppler signal in color flow mapping: • The frequency characteristic of the estimator, which indicates the range of frequencies for which the relation between the estimated and the true mean frequency of the Doppler signal are reasonably equal, and which is directly linked to the range of flow velocities that the device can manage. • The robustness of the mean frequency estimation
905
in the presence of noise, which directly determines the quality of Doppler images. These parameters are studied for an ideal noiseless signal and for a signal-to-noise ratio ( S / N ) of 0 dB. These values have been chosen from the results of a more extensive study. They are representative of the behavior of estimators when analyzing excellent and rather poor Doppler signals.
Doppler signal model The Doppler signal has been modeled by a complex value signal whose real and imaginary parts are in quadrature to comply with the requirements of the baseband conversion. The Doppler spectrum is rectangular, and its bandwidth can be continuously set between AB = 0 and AB = PRF. This deterministic model of the Doppler signal has been preferred to a stochastic model. It will allow to separate clearly, in the subsequent studies, the respective contribution of the Doppler signal and that of the noise. However, a complementary study, not reported here, and using the Doppler signal model proposed by Mo and Cobbold (1986), has shown that the conclusions provided in this article are comparable for both models. The nature of the noise physically existing on the Doppler signal is hard to model. The real and imaginary components of the electronic noise are mutually uncorrelated. The baseband conversion of the radiofrequency signal leads to a phase difference of 7r/2 between the two components of the noise. Finally, cross-talk phenomena contribute to correlation between these two components. The behavior of the estimators as a function of the noise has thus been studied for white Gaussian and uniformly distributed noises. For each type of noise, equal, in quadrature or completely uncorrelated, real and imaginary parts have been used, and the results obtained in all cases remain identical. The present study corresponds to a white, evenly distributed noise with uncorrelated real and imaginary parts. The frequency characteristics of estimators (mean frequency estimated as a function of the real mean frequency of the Doppler signal) have been obtained by shifting the mean frequency of this signal between - P R F / 2 and + P R F / 2 . The comparison between the behavior of both estimators is achieved for an N = 48 Doppler sample window to eliminate the influence of the signal windowing on the estimation of correlation lags and the truncation effect in the sum of expression (24).
Frequency characteristic Fig. 2 displays the frequency characteristic of
CE~ and that of CEn for a Doppler signal bandwidth
906
Ultrasound in Medicine and Biology
+PRF/2
~'
-PRF/2
Volume 17, Number 9, 1991
+PRF/2
a
i
I
I
A
..o.°o'°
I
+PRF/2
+PRF/2
-PRF/2
•
i
i
i
B
+PRF/2
D
+PRF/2
+PRF/2
°,~oqO°~..~*o°°°
~o ,.°°°°° -PRF/2
C
+PRF/2
-PRF/2
Fig. 2. Frequency characteristics of C E I and CEn:x = actual Doppler frequency, y = estimated Doppler frequency. CEI:dotted line, CE,: plain line. (A) N = 48, AB = 0.001 PRF; (B) N = 48, AB = 0.2 PRF; (C) N = 8, AB = 0.001 PRF; (D) N = 8, AB = 0.2 PRF. Note the inflexions of the frequency characteristic corresponding to _+(PRF-AB)/2, here _+0.4 PRF for CEn in Fig. 2(B).
AB = 0.001 PRF (very narrow band signal) and AB = 0.2 PRF. It clearly shows that CE~ has a frequency range covering the whole Nyquist domain and that CE, has a frequency range that directly depends on the bandwidth of the Doppler signal. These results are in accordance with the theoretical behavior of both estimators (Kristoffersen 1988): The assumption of a continuous time signal does not a priori impose constraints on the frequency range of CEI, which will be able to estimate mean frequencies all along the Nyquist band from - P R F / 2 to + P R F / 2 as far as the phase ~ 1) does not vary too rapidly. In the present case, this estimate is unbiased because the rectangular spectrum chosen is symmetrical with respect to its mean frequency. The assumption of a discrete time Doppler signal and the use of a development in Fourier series of a periodic PSD imply that the mean frequency estimated by CEn will only be significant when the lowest
frequency of the Doppler spectrum will be higher than - P R F / 2 and when its highest frequency will be lower than +PRF/2. Therefore, the frequency range of CE, will be equal to [PRF - AB]. Outside these limits, the bias increases very rapidly along with frequency. Figure 3 presents the average of 64 of the previous frequency characteristics for S / N = 0 dB and AB = 0.2 PRF. One can observe a reduction in the frequency range for both estimators. This can be explained by the increase in the estimation variance when the extreme frequencies of the spectrum cross _+PRF/2. As a result, the averaging of mean frequency, which is commonly used in color flow mapping to decrease the noise on the image, is not really adequate for noisy signals. Indeed, it induces a decrease in the frequency range, which directly depends on the spectral width of the Doppler signal and on the signal-to-noise ratio.
Adaptive estimation for Doppler signal • A. HERMENTel al.
+PRF/2
.
907
+PRF/2
.
.
,
.
-PRF/2
,
•
A
.
+PRF/2
+PRF/2
:
:
;
-PRF/2
.
.
.
.
B
+PRF/2
+PRF/2
.
-PRF/2
:
C
.
.
.
I
+PRF/2
-PRF/2
I
I
.
I
D
.
.
.
+PRF/2
Fig. 3. Average m e a n frequency characteristic of CE, and CEn for a noisy Doppler signal, S / N = 0 dB, AB = 0.2 PRF: x = actual Doppler frequency, y = averaged estimated Doppler frequency. ( A ) CE,, N = 48; (B) CE~, N = 48; ( C ) CE~, N = 8; ( D ) CE,, N = 8.
Estimation variance Figure 4 presents the standard deviation (SD) of a mean frequency estimated as a function of the signal mean frequency, for AB = 0.2 PRF and S / N = 0 dB. Figure 5 presents the SD for these two estimators as a function of the signal bandwidth for S / N = 0 dB, but for a fixed zero mean frequency. Again, these results are in accordance with the theoretical behavior of both estimators (Kristoffersen 1988), in particular, the rather constant variance of CE1 as far as no aliasing occurs, and the strong increase in this variance when the extreme frequencies of the Doppler spectrum cross the Nyquist limits. CEn has a far better noise immunity than CEI, and this is even more observed as the signal bandwidth increases. The performances of the estimators directly depend on their ability to estimate the derivative at the origin of the correlation function. If the correlation function varies slowly around the origin, the correla-
tion lags evolve slowly. The derivative is correctly estimated from C( 1 ) alone for CEI and from the sum of eqn (24) for CErt. If the variations of the correlation function at the origin become faster (high mean frequency or wideband signal), the derivative will be less and less well estimated by CEI. On the other hand, the correlation lags decrease faster, the sum ofeqn (24) will converge faster, and CEn will adequately estimate the mean frequency in the [ - ( P R F - AB)/2, + ( P R F - AB)/2] range, thus accounting for a lower variance. In addition, eqn (14) shows that CE~ only uses the first correlation lag, and any uncertainty on the value of Cz( 1 ) is directly converted into an uncertainty on the mean frequency. The impact of this uncertainty will be minimal when Im{ Cz( 1 )} will be small and when Re{ Cz(1)} will be large, i.e., to say for the low mean frequencies; conversely, it will be maximal around __+PRF/2. On the other hand, CEn utilizes a weighted aver-
908
Ultrasound in Medicine and Biology +PRF/2
Volume17, Number 9, 1991 +PRF/2
'
lr
i
-PRF/2
A
+PRF/2
+PRF/2
:
'
I
e
I
I
B
I
+PRF/2
+PRF/2
., -PRF/2
'
;
-PRF/2
C
+PRF/2
-PRF/2
:
_,
:
:
D
;
:
:
+PRF/2
Fig. 4. Standard deviation of CEI and CE, as a function of the signal mean frequency, S/N = 0 dB, AB = 0.2: x = actual Doppler frequency, y = standard deviation of the estimated Doppler frequency. (A) CE~, N = 48; (B) CE,, N = 48; (C) CE~, N = 8; (D) CE,, N = 8.
age of the Im{ Cz(n) }, and the uncertainties in calculation o f the Cz (n) may be expected to decrease in the summation of eqn (24).
Influence of the truncation of the Doppler signal T o comply with real-time color Doppler imaging requirements, the n u m b e r of samples available to calculate the mean frequency is reduced ( 4 - 1 2 samples). In the following, if not otherwise mentioned, Nwill be taken equal to 8 for short windows. This will have a direct impact on the sum in eqn (24) for CE, and on the estimation of correlation lags for both estimators.
Influence of truncation on the mean frequency expression The truncation on the sum in eqn (24) will lead to an increased Gibbs phenomenon, which will generate oscillations in the frequency characteristic, even if the correlation lags are perfectly estimated. This phenomenon is illustrated in Fig. 6. It displays the fre-
quency characteristic of CE, for the narrow band signal (zkB = 0.001 ) and when an N = 256 window is used for correlation lags computation. Figures 6 (A,B) corresponds to the utilization o f the only first 48 and the only first 8 of these previous correlation lags in eqn (24), respectively. The oscillations in the frequency characteristic are even more important as the number of lags is low.
Influence of the window length on the frequency characteristic Figure 2 allows one to compare the frequency characteristics of CE, and that of CEn for long and short windows for 2xB -- 0.001 PRF and a d / = 0.2 PRF. Comparison of these curves shows that the frequency characteristic of CE1 is independent of the window length and that the one o f CE, depends on the window length, besides the already observed dependence of this frequency characteristic with the signal bandwidth for CE,.
Adaptive estimation for Doppler signal • A. HERMENTet aL
909
0,3
// 0,2
0,1
B'----~
0,0
0,0
r-
0,2
0,4
0,6
0,8
1,0
Fig. 5. Standard deviation of the mean frequency estimation by CE~ and CE, for a noisy signal as a function of spectrum width; S/N = 0 dB: x = actual Doppler bandwidth, y = standard deviation of the estimated Doppler frequency. (A) CE,, N = 48; (B) CE,, N = 48; (C) CE~, N = 8; (D) CE,, N = 8.
Figure 3 compares the average m e a n frequency characteristics of CE~ and CE, for long and short windows, for AB = 0.2 P R F and for a 0 dB signal-to-noise ratio. The general influence of the decrease in the window length is a m o r e i m p o r t a n t degradation of the average frequency characteristic for CE~ than for CE,, even though the frequency range o f CE~ is always better than the one o f CE, for a given window.
+PRF/2
Influence of the window length on the estimation variance Figure 4 indicates the modification o f the SD o f m e a n frequency estimation o f CE~ and CE, along with the frequency, for long and short windows, a n d for S / N = 0 dB. For both estimators, the SD increases with the reduction of the window length. For longer windows, the correlation function o f
+PRF/2
f
/
-PRF/2
A
+PRF/2
-PRF/2
B
+PRF/2
Fig. 6. Influence of a truncation in eq. (24), N = 256, AB = 0.001 PRF: x = actual Doppler frequency, y = estimated Doppler frequency. (A) 48 correlation lags; (B) 8 correlation lags.
910
Ultrasound in Medicine and Biology
the noisy part of the signal tends toward the theoretical correlation function of the noise, which has no contribution on the correlation lags excepted on Cz(0). In addition, the SD of CEn remains rather constant with frequency, whatever the window when that of CE~ is degraded around + / - PRF/2 for short windows as a result of the CE~ behavior itself and of an artificial enlargement of AB caused by window truncation. Figure 5 compares the SD of these estimators for the same 0 dB signal-to-noise ratio and the same N = 8 window. Again, one can observe that the SD increases with the decrease of N, CE~ being more sensitive to the window reduction. Moreover, the SD remains independent of the signal bandwidth for CEn, while it is extremely AB-dependent, tending toward a constant value when the signal bandwidth increases, for CE1, whatever the N value. Starting from the results of this comparison, the question of defining an estimator allowing to conciliate a high-frequency range and a low-estimation variance is now raised. For that purpose, the starting point will be CEn given by eqn (24), CE~ being completely defined by eqn (13).
Optimization of CEn Three solutions have been studied to improve the
CEn estimator: first, to replace the weighting coefficients an = (-l)"+~/a'n by coefficients taking into account the Gibbs phenomenon associated with the expansion truncation; second, to extrapolate the correlation function so as to create a greater number of correlation lags; and, third, to interpolate the correlation function for the same purpose.
Modification of a, coefficients Here, the problem consists in finding the best polynomial approximation of the ideal frequency characteristic. It has been shown that using the coefficients an = (-1 )n + l / r n associated with the unbiased correlation estimator C2 of eqn (30) led to a larger frequency range but to a more oscillatin~ characteristic of the estimator. In fact, changing C~ to C2 gives the same estimator equation as changing the theoretical weighting sequence { an } n = 1, 2 , . . . , N - 1 for the following:
Volume 17, Number 9, 1991
prove the frequency range and increase the Gibbs phenomenon.
Exponential weighting A simple means to reinforce the frequency characteristic of the estimator is to choose a weighting sequence such as: ( - 1 ) n+,
an -
an
( - 1 ) n+l
N- n
lrn
(32)
with ~ < 1. Figure 7 displays, for N = 8, for AB = 0.001 PRF, and AB = 0.4 PRF, a comparison between the estimator CEn of eqn (24) and the exponentially weighted estimator defined in eqn (32) for/3 = 0.6. The absolute maximum value of the estimated mean frequency increases along with the decrease of /3, but the frequency range remains unchanged. The frequency range of the new estimator remains, in fact, directly connected to the bandpass of the Doppler signal as for CEn.
Least squares weighting A more subtle approach has been proposed by Angelsen and Kristoffersen (1983a). It consists in determining the series of coefficients dn so as to minimize the mean quadratic error on:
J,=o that is to say, to calculate the dn coefficients minimizing:
LI L d.,- J,.=o L d,-J,-oolJ
(33)
where E represents the expectation operator, Cz(z) the estimated correlation function and Cz(r) the true correlation function. For that purpose, one must define a probability law for the PSD obtained in order to calculate the previous mean quadratic error:
---&---, j,=o - t---&-, J,=o =
1
71-na
Dk(f)gz(f)df
(34)
(31)
while keeping the same C~ estimator. One can argue that any change in the sequence { an } that will increase the relative weight values of an with n, will i m -
where G(f) represents the frequency band of the Doppler signal and where: k
Dk(f) = X [d,e -z'f'] - 2iTrf n=-k
(35)
Adaptive estimation for Doppler signal • A. HERMENTet al.
+PRF/2
.
911
+PRF/2
.
.
i
.
-PRF/2
;
A
;
I
+PRF/2
+PRF/2
i
i
I
i
-PRF/2
I
I
B
I
+PRFI2
+PRF/2
I
-PRFI2
C
+PRF/2
-PRF/2
I
I
D
I
+PRF/2
Fig. 7. Influence of an exponential weighting of the a, series on CEn, N = 8. x = actual Doppler frequency, y = estimated Doppler frequency. (A) A B = 0.001 PRF;/~ = I;(B) A B = 0.001 PRF, ~ = 0.6;(C) A B = 0.4 PRF;/~ = l; (D) A B = 0.4 PRF, B = 0.6.
By choosing an even probability law on G(f) and a line spectrum gz(f) = ~( f ) , the minimizing criterion is then written as follows:
min
d,, = fc
IDk(f)i2df
+PRF/2
A
,J"
(36)
•/
\
(f)
Therefore, the dn are the coefficients o f the trigonometric polynomial that is closest, to 2 h r f o n G ( f ) , in the least squares sense. One can note that if G(f) = [ - P R F / 2 , + P R F / 2 ] , one again finds the coefficients an = (-1)"+1/2mr, which are the Euler coefficients for expanding 2i~-fin a Taylor series. In this approach, if G(f) reduces to [ - f 0 , + f 0 ] , the approximation will be better on G(f) iff0 is low for a given n. Figure 8 displays a comparison, for AB = 0.001 PRF, N = 8 samples and for G(f) = [ - 0 . 4 5 PRF, +0.45 P R F ] , between the frequency characteristics o f CE, in dotted lines and the one corresponding to the
j,"
f
j
"~
t"' '/ ."
-PRF/2
~
I
|
0
I
!
I
+PRF/2
Fig. 8. Influence of a least squares optimization of the an series, N = 8, AB = 0.001 PRF: x = actual Doppler frequency, y = estimated Doppler frequency. (A) an = ( - 1) n+, / ~rn; (B) an coefficients according to Angelsen and Kristoffersen (1983a).
912
Ultrasound in Medicine and Biology
coefficients given by Angelsen and Kristoffersen (1983a), as a broken line. The optimized characteristic of the estimator is more linear, but its frequency range is still unsatisfactory. This approach does not either seem to Satisfy the requirements as originally stated. Kristoffersen (1988) applied the Same least squares approach to design a mean frequency estimator offering a constant ~f/ferror characteristics. The result of this approach is comparable to the previous exponential weighting: the maximum value of the curve increases, but its frequency abscissa remains the same.
Extrapolation of the correlation function for n > N Signal modeling allows one to extrapolate the correlation function for n > N. The most direct method (Kay and Marple 1981 ) models the correlation function Cz(n) by an AR model of adequate order, and extrapolates Cz(n) by a Burg algorithm that allows one to define, among all the possible extrapolations of the correlation function, the one that has a maximum entropy (Schwindlein and Evans 1989). Beyond N, however, extrapolation must be limited because the stability of estimation cannot always be ensured. A good compromise would consist of taking, for Doppler signal analysis, [0 < n < 2N] (Schwindlein and Evans 1989). However, the complexity of this approach does not bring a significant gain on the estimation of mean frequency. The extrapolated values of the correlation
function for n > N have a weak ( 1/~rn) contribution in eqn ( 24 ). Last, the order (P) of AR models must be limited to P < N for an estimation of this model by least squares techniques. For Doppler signals, it is generally impossible to obtain a ratio higher than P / N = 1/4 (Kaluzinski 1989). One must notice that color flow mapping leads to 4 > N > 12, i.e., AR model orders smaller than 4, which do not allow for a correct Doppler signal model. This model order limitation can be removed with long AR models (P = N) associated with regularized least squares methods (Demoment et al. 1988). The regularization principle, applied here, consisted in estimating the AR parameters while minimizing a distance that not only depends on the fidelity of the model to the data (least squares approach), but also on a quadratic criterion characterizing the AR spectrum smoothness. The extrapolated series of correlation lags is then obtained from AR coefficients using an inverse Levinson algorithm. Despite this additional effort, the results obtained are still almost equivalent to those directly provided by eqn (24). The frequency range is still limited by the window length as shown in Fig. 9 for AB = 0.001 PRF and AB = 0.4 PRF. In the present case, it is likely that the increase in the model order contributes in itself to a better description of the Doppler signal but that, on the other hand, the AR coefficients obtained with regularization techniques are mutually correlated and that the overall result remains equivalent to that of the previous approach. To conclude, the complexity introduced by these methods does not seem to be justified with respect to +PRF/2
+PRF/2
-PRF/2
Volume 17, Number 9, 1991
A
+PRF/2
-PRF/2
B
+PRF/2
Fig. 9. Correlation lag estimation by AR modelling, N = 8: x = actual Doppler frequency, y = estimated Doppler frequency. (A) AB = 0.001 PRF; (B) AB = 0.4 PRF.
Adaptive estimation for Doppler signal Q A. HERMENT el al.
913
the possible improvement it would bring to the mean frequency estimation.
ues C~(n) = C~(3m) by writing Cz(n) again in the following form:
Interpolation of correlation lags An interpolation of the correlation lags C~(n) will not only lead to more available correlation lags to estimate f, but also to create new intermediate correlation lags corresponding to a higher sampling frequency, because their delay is a submultiple of the sampling period. This will not eliminate aliasing problems beyond PRF/2, but it will enable the operation on the linear part of the frequency characteristic. For instance, creating a correlation lag interpolated between each correlation lag Cz(n) would lead to an equivalent N = 16 window and only use the f values from - P R F / 4 to + P R F / 4 on the frequency characteristic. However, as the vertical scale is doubled because of oversampling, this will correspond to a real range of - P R F / 2 to + P R F / 2 while keeping the linearity of the center part of the frequency characteristic of CE, between - P R F / 4 and + P R F / 4 . The first option consists in defining the number of correlation lags to interpolate between successive original correlation lags C, (n). Trials have been made on different orders. It appears that the addition of two lags interpolated between two consecutive Cz(n) allows, for short windows, optimization of the ratio between the quality of estimation of the mean frequency and the amount of additional calculations introduced by interpolation. The second option consists in determining the interpolation scheme. To obtain an estimator with a high-frequency range, we have chosen to elaborate a scheme allowing to get closer to hypotheses underlying to the setting up of estimator CE~, i.e., an interpolation that tends to reduce the variation of the modulus of the correlation function at the origin and that strengthens the influence to the phase of the first correlation lag. For that purpose, the modulus and the phase of the interpolated correlation lag were linearly interpolated between the successive original correlation lags. One would note that this interpolation is fundamentally different from a linear interpolation of the real and imaginary parts of the Cz(n). The interpolated estimator CE3n is thus defined by:
Cz(n) = ICz(n)le i~")
2
f = l -- N-| E [(-l)3m-t+la3m-tlm{C'~(3 m - l ) } ] 71" m=l /=0
(37) The values ofC'( 3m - l) are calculated from the val-
(38)
defined by its modulus ICz(n) I and its argument ~o(n). The interpolation scheme is: IC'z(3m)l = ICz(n)l
(39)
I C ~ ( 3 m - 1)1 = [2[Cz(n)l + I C z ( n - 1)1]/3
(40)
=[[Cz(n)l+2lCz(n-1)l]/3
(41)
C'z(3m - 2)1
~0'(3m) = ~o(n)
(42)
~'(3m - 1 ) = [2~o(n) + ~o(n - 1)1/3
(43)
~o'(3m - 2) = [~0(n) + 2~o(n - 1)]/3
(44)
C'z(3m) = IC'~(3m)le i¢~3m)
(45)
C'z(3m- 1 ) = IC'z(3m- 1)l e~¢t3'~-l)
(46)
C'z(3m - 2) = IC~(3m - 2)le i¢~3"-2)
(47)
The a3m_l coefficients are also interpolated by: aam = 3 / 3 m 1)
(49)
a3m_2 = 3 / ( 3 m - 2)
(50)
a3m_ 1 =
3/(3m-
(48)
The frequency characteristic of CE3n is identical to that of CE~ for AB < 0.4 and slightly differs for wider spectra. However, as illustrated in Fig. 10, even for large bandwidth signals (AB = 0.9 PRF) and a short window, there is an obvious similarity between the frequency characteristics of both estimators. As for the estimation variance, it is strictly equivalent for both CE3n and CE~. Finally, note that ifa3m_l = a3m-2 = 0 , then CE3n = CE~. Hybrid estimator Equations (37-50) allow for a transition from estimator CEn to estimator CE1. Therefore, it seems possible to define, by an appropriate weighting of coefficients a3m-1 of eqns (48-50), an estimator halfway between CEn and CE~ and, in particular, to adapt the compromise between the frequency range and the variance of estimation to the characteristics of the sig-
914
Ultrasound in Medicine and Biology
+PRF/2
+PRF/2
-PRF/2
Fig. 10. Comparison between CE, and C&, N = 8, AB = 0.9 PRF: x = actual Doppler frequency, y = estimated Doppler frequency. Dotted line: CE, ; plain line: CE,,
nal to process. We have thus generated a hybrid estimator CE,,, defined as follows:
f=
b2;[(-1)3”-‘+‘aJ,_,(cX)
Volume 17, Number 9, 199 I
However, if some characteristics on the Doppler equipment, such as the noise level of the system are known a priori, or if some information on ‘the Doppler signal characteristics can be inferred from the environment of the considered window, the characteristics of the mean frequency estimator can be adapted, for each window, to process the corresponding signal. The following scheme presents one simple application of this adaptive estimator to the Doppler signal. The adaptivity is based on the assumption that the Doppler signal mean frequency cannot be very different between two adjacent windows:
.L+,%fk
(55)
In the case of a monogate or TM Doppler system, k will characterize two windows located at the same site but at consecutive instants. For a color Doppler imaging system, k will correspond to two windows on different sites: either two consecutive windows along the same scanning line, or two windows located at the same depth along two consecutive scanning lines. The criterion selected to simply express the gradual changes of the Doppler signal is as follows:
I=0 X Im{C;(3m
-
f)j]
(51)
The values of C:( 3m - I) are calculated with eqns (38,47), but a3m_1coefficients are replaced by new coefficients u3,,_,( LY)defined as follows: ~~~(a) = (1 +
2&)/3m
(52)
a3m-1(cu) = a( 1 + 2a)/( 3m - 1)
(53)
= (Y(1 + 2a)/(3m - 2)
(54)
a3&~)
(Yis a trimming parameter ranging between [ 0 and l] and determined by the user.
Adaptive estimator An improvement can thus be obtained for a given signal by choosing the value (Ythat is most appropriate to the characteristics of this signal. Unfortunately, the statistical properties of the Doppler signal change very rapidly throughout the cardiac cycle, even for a given location of the flow. The signal-tonoise ratio strongly varies between diastole and systole as a result of hemodynamical phenomena inherent in flow conditions and to the echo canceller cutoff characteristics. Moreover, the signal bandwidth can also be strongly modified according to flow conditions during the same heart contraction, which leads to a difficult choice for (Y.
(56) a2 allows one to smooth the variation of fbetween Sk and Jik+,. A low value for a2 impedes any rapid modification of the mean frequency between two consecutive signal windows. As for the coefficient QI, , it avoids an oversized CJ (Y, , a2) in eqn (56) when the mean frequency is close to zero and permits mean frequency variations when fk = 0. An illustration of the respective contribution of these two parameters is given in Fig. 11. The detailed block diagram of Fig. 12 displays an algorithm that chooses the largest (Yvalue complying with the smoothness criterion Cr( a,, a2). a0 represents the incremental step of CL
al
t&small
a2 too small
Fig. 11. Influence of a, and (Yeon the mean frequency estimation. x = time; y = estimated Doppler frequency.
Adaptiveestimationfor Dopplersignal• A. HERMENTet al.
915
+PRF/2
k=0 f0=0 OCO,(gl, 0[2
a=l+ao
k=k+l
I
k-
•
I
l
I
•
•
•
+PRF/2
Fig. 14. Average mean frequency characteristic of CE3,,,,, CE~ and CE,, for a noisy Doppler signal, S / N = 0 dB, AB = 0.2 PRF, N = 8: x = actual Doppler frequency, y = averaged estimated Doppler frequency. (A) CEI; (B)
f = EC3n~
fEn;
F,
I
-PRF/2
(C) CE3na,
ot =
0.8.
Cf ( a l , tI2)
lno
~=0
RESULTS
I no
~ i Y es k=k~
I no
yes
e
END
I
Fig. 12. Principle of adaptive mean frequency estimation with CE3n~.
Properties of the hybrid estimator Figures 13-16 display the b e h a v i o r o f CE3n . for a = 0.8, N = 8 a n d AB = 0.2 P R F together with those o f CE~ a n d CE, for c o m p a r i s o n o f their respective performances. Analysis o f these curves does show, as expected, for CE3n,~, an intermediate b e h a v i o r between that o f CE t a n d that o f CE, b o t h for the bias a n d the variance. A n y a value larger t h a n 0.8 will give to CE3,,, a behavior closer to that o f CEI a n d a n y a value lower t h a n 0.8 a behavior closer to that o f CE,.
+PRF/2 +PRF/2
I
-PRF/2
+PRF/2
Fig. 13. Frequency characteristics of CE3n,,, CEt and CE,,, N = 8, AB = 0.2; PRF: x = actual Doppler frequency, y = estimated Doppler frequency. (A) CE~; (B) CE,,; (C) CE3.,,, a = 0.8.
-PRF/2
I
I
I
+PRF/2
Fig. 15. Standard deviation of C3,~, CEi and CE, as a function of the signal mean frequency, S / N = 0 dB, AB = 0.2 PRF, N = 8: x = actual Doppler frequency, y = standard deviation of the estimated Doppler frequency. (A) CE~ ; (B) CE.; (C) CE3,,,,, a = 0.8.
916
Ultrasound in Medicine and Biology
0,2
Volume 17, Number 9, 1991
0,6
094
098
190
Fig. 16. Standard deviation of the mean frequency estimation of CE,,,. CE, and CE, for a noisy signal as a function of spectrum width; N = 8, S/N = 0 dB: x = actual Doppler bandwidth, y = standard deviation of the estimated Doppler frequency. (A) CE, ; (B) CE,,; (C) CE,,,, (Y= 0.8.
Application of the adaptive estimator to actual Doppler signals Doppler signal recordings have been achieved using a pulsed velocimeter operating at 3 MHz. Two sample volumes were located downstream from a stenosis, one in the laminar jet and the other in the area of turbulent flow corresponding to the breakup zone
of the jet. These recordings have been made in vitro using a pulsatile flow. In the following application, N = 4 sample windows were used to test the adaptive estimator under stringent conditions. The gap between two consecutive windows k and k + 1 was of 60 sampling periods as explained in Fig. 17.
1--11111-111-11
N=60
N=4
Doppler window Tk
N=60
Unprocessed data
N=4
N=60
Doppler wicidow fk+l
Fig. 17. Timing used for adaptive processing of actual Doppler signals. x = time, y = Doppler signal amplitude.
Adaptive estimation for Doppler signal • A. HERMENTel al.
+PRF/2
917
+PRF/2
1 •"
-PRF/2
|
l
A
0.75 s
+PRF/2
-PRF/2
B
+0.75 s
+PRF/2
'I
-PRF/2
C
0.75 s
t-PRF/2
|
0/ .1
ii
|,,,
D
0.75 s
+PRF/2
Fig. 18. Behavior of CE~, CE, and CE3,,, on Doppler signals: laminar flow. x = actual Doppler frequency, y = estimated Doppler frequency. (A) CE~, N = 48; (B) CEt, N = 4; (C) CE,,, N = 4; (D) CEa,,,,, N = 4; (E) CE,,, premixing, N = 4.
-PRF/2
E
0.75 s
The results, obtained using the algorithm described in Fig. 12, are compared with those of C E i and CE,, for laminar flows on Fig. 18 and for a turbulent flow on Fig. 19. On these two figures, curves A correspond to the mean frequency estimation using CE~ and an N = 48 sample window. They will be considered as the exact mean frequency of the Doppler signal even if this is a rough approximation
in the case of a turbulent flow. Curves B, C, D and E correspond to CE~ , CE,, CE3n a and CE,,, when using the premixing scheme proposed by Kristoffersen (1985). In use of the adaptive estimator, the same values at = 0.001 and a2 = 0.5 were used for both laminar and turbulent flows. In this premixing scheme, the frequency shift of the Doppler signal was simply taken
918
Ultrasound in Medicine and Biology +PRF/2
Volume 17, Number 9, 1991 +PRF/2
r-
""]
-PRF/2
A
r
0,75 s
I'!t -PRF/2
I' lpl B
+0.75 s
D
0.75 s
+PRF/2
+PRF/2
-PRF/2
C
0.75 s
-PRF/2
+3 PRF/2
Fig. 19. Behavior of CEI, CE, and CE3,~ on Doppler signals: turbulent flow. x = actual Doppler frequency, y = estimated Doppler frequency. (A) CE~, N = 48; (B) CE,, N = 4; (C) CE,,, N = 4; (D) CE3,,a,N = 4. (E) CE,,, premixing, N = 4.
O...l
~r'-~ -PRF/2
I E
0.75 s
equal to the mean frequency fk. estimated in the previous time window. The adaptivity of the CE3n,~estimator allows one to improve significantly the quality of estimation of the mean frequency by adequately estimating the highest frequencies of the systolic phase while keeping a reduced estimation variance during the diastole.
When the signal premixing down to the zero frequency is applied to the laminar flow Doppler signal, the mean frequency is excellent. However, for the turbulent flow, an important mean frequency error propagates all along the cycle because o f the transmission of local erroneous mean frequency estimations. From these results, improvement of color flow
Adaptive estimation for Doppler signal • A. HERMENTel aL
mapping images can be contemplated. Adaptive mean frequency estimation should improve the image quality for a given frame rate, but also increase the frame rate while keeping comparable images. This should be of some interest for cardiac imaging where the flow modifications are extremely fast during the systole. Finally, it can offer a partial contribution to improvement in the compromise between low-frequency cutoff and rise time in the echo cancellers. The excellent noise immunity of the adaptive estimator for low frequencies, as for the CE,,, should allow better analysis of low frequencies. CONCLUSION The color Doppler mean frequency estimator
CE~ and the general mean frequency estimator CE, are based on the evaluation of the first derivative of the Doppler signal correlation function at the origin. The first one has a frequency range of [ - P R F / 2 , + P R F / 2 ], but its variance of estimation is important and its noise immunity very poor. The second one has a reduced frequency range but a low estimation variance and an attractive noise immunity. Using a nonlinear interpolation of the different correlation lags for the Doppler signal correlation function derivative estimation, has led to a new mean frequency estimator. Its behavior is continuously shifted from that of the CE~to that of the CE,, according to the value of a single trimming parameter. This hybrid estimator allows one to optimize, for a given window, the compromise between the frequency range and the estimation variance to the particular features of a given Doppler signal. Its preliminary adaptive application to the Doppler signals reveals its interest when one is processing very short Doppler windows. However, the best adaptivity criterion has still to be defined for each particular application, in particular, adaptivity schemes fitted to color Doppler imaging. Acknowledgments--We gratefully acknowledge Mrs. F. Marchand for her editorial and secretarial assistance.
REFERENCES Angelsen, B. A.; Kristoffersen, K. Discrete time estimation of the mean Doppler frequency in ultrasonic blood velocity measurements. IEEE Trans. Biomed. Engin. BME-30:207-214; 1983a. Angelsen, B. A.; Kristoffersen, K. Zero-crossing density for ultrasonic Doppler signals obtained from computer simulations. Ultrasound Med Biol 9:661-665; 1983b.
919
Arts, A. G. J.; Roevros, J. M. On the instantaneous measurement of blood flow by ultrasonic means. Med. Biol. Eng. 10:23-34; 1972. Barber, W. D.; Eberhard, J. W.; Carr, S. G. A new time-domain technique for velocity measurements using Doppler ultrasound. IEEE Trans. Biomed. Eng. BME-32:213-229; 1985. Bramanti, M.; Marchesini, E. Analysis ofa PLL system as analog processor in CW Doppler flowmeters. IEEE Trans. Biomed. Eng. BME-3:584-589; 1983. Demoment, G.; Houacine, A.; Herment, A.; Mouttapa, I. Adaptive Bayesian spectrum estimation. Fourth ASSP Workshop on spectrum estimation and modeling, Minneapolis, USA, August 3-5; 1988. Gerzberg, L.; Meindl, J. D. Power spectrum centroid detection for Doppler systems applications. Ultrason. Imaging 2:232-261; 1980a. Gerzberg, L.; Meindl, J. D. The ~ power spectrum centroid detector: System considerations, implementation and performance. Ultrason. Imaging 2:262-289; 1980b. Kaluzinski, K. Order selection in Doppler blood flowsignal spectral analysis using autoregressive modelling. Med. Biol. Eng. Comput. 27:89-92; 1989. Kasai, C.; Namekawa, K.; Kouano, A.; Omoto, R. Real-time twodimensional blood tlow imaging using an autocorrelation technique. IEEE Trans. Sorties Ultrasonics SU-32:458-464; 1985. Kay, S.; Marple, S. L. Spectrum analysis. A modern perspective. Proc. IEEE 69:1380-1419; 1981. Kristoffersen, K. Time domain estimation of the center frequency and spread of the Doppler spectra in diagnostic ultrasound. IEEE Trans. Ultrason. Ferroelec. Freq. Contr 35:685-700; 1988. Kristoffersen, K.; Angelsen, B. A. J. A comparison between mean frequency estimators for multigated Doppler systems with serial signal processing. IEEE Trans. Biomed. Eng. BME 32:645-657; 1985. Kuc, R.; Li, H. Reduced-order autoregressive modeling for centerfrequency estimation. Ultrason. Imaging 7:244-251; 1985. Loupas, T.; McDicken, N. Low order complex AR models for mean and maximum frequency estimation in the context of Doppler color flow mapping. IEEE Trans. Ultrason. Ferroelec. Freq. Contr 37:590--601; 1980. Miller, K. S.; Rochwarger, M. M. A covariance approach to spectral moment estimation. IEEE Trans. Information Theory 18:588596; 1972. Mo, L. Y.; Cobbold, R. S. Speckle in continuous wave Doppler ultrasound spectra: A simulation study. IEEE Trans. Ultrason. Ferroelec. Freq. Cont. 33:747-753; 1986. Nowicki, A.; Reid, J. M. An infinite gate pulsed Doppler. Ultrasound Med. Biol. 7:41-50; 1981. Reid, J. M. Processing and display techniques for Doppler flow signals. In: Harrison, D. C.; Sandier, H.; Miller, H. A., eds. Cardiovascular image and image processing. The society of photooptical instrumentationEngineers; 1966:73-78. Rice, S. O. Mathematical analysis of random noise. BSTJ 23:288; 1963. Sainz, A.; Roberts, V. C.; Pinardi, G. Phase-locked loop techniques applied to ultrasound Doppler signal processing. Ultrasonics 14:128-132; 1976. Schlindwein, F. S.; Evans, D. H. A real-time autoregressive spectrum analysis for Doppler ultrasound signals. Ultrasound Med. Biol. 15:263-272; 1989. Talhami, H. E.; Kitney, R. I. Maximum likelihood frequency tracking of the audio pulsed Doppler ultrasound using a Kalman filter. Ultrasound Med. Biol. 14:599-609; 1988.
0301-5629/91 $3.00 + .00 © 1991 Pergamon Press plc
Printed in the U.S.A.
OOriginal Contribution ADAPTIVE ESTIMATION OF THE M E A N FREQUENCY OF A DOPPLER SIGNAL FROM SHORT DATA W I N D O W S A. HERMENT, t G. DEMOMENT,* J. P. GUGLIELMI, t PH. DUMEE t a n d C. PELLOTt tINSERM U.256, H6pital Broussais, Paris, France; *L2S(CNRS/ESE), UPS, Gif-sur-Yvette, France (Received 7 March 1991; in final form 8 August 1991 )
Abstract--The color Doppler estimator ( C E 1), which is calculated from the phase of the first correlation lag of the Doppler signal, is compared to the general mean frequency estimator (CE.), which is based on a weighted summation of all the available correlation lags, for long and short Doppler data sets (typically 48 and 8 Doppler samples). A new estimator of the Doppler signal mean frequency is derived from the results of this study. It optimizes the compromise between the range of analyzable frequencies and the estimation variance for the characteristics of the Doppler signal. Demonstration is provided that the behavior of this estimator shifts from that of CEI to that of CE,, according to the setting of a single parameter. An adaptive version of this estimator is implemented and applied to Doppler recordings. Applications can be contemplated for color Doppler imaging.
Key Words: Color Doppler, Mean frequency, Correlation function.
INTRODUCTION
ties of the correlation function of the Doppler signal, in particular, the double correlator (Gerzberg and Meindl 1980a) and the single correlator (Arts and Roevros 1972). This article presents a comparison between two of these correlation-based estimators of the Doppler signal mean frequency. The first estimator is the conventional color Doppler estimator (CE~), introduced by Miller and Rochwarger (1972), and also referred to as the correlation angle estimator, which utilizes the phase of the first correlation lag of the Doppler signal. The second estimator is the general mean frequency estimator (CE,), described by Angelsen and Kristoffersen (1983a), and also referred to as the mean frequency estimator, which utilizes a weighted sum of all the available correlation lags of the Doppler signal. This study will be used as a starting point to new adaptive mean frequency estimator of the Doppler signal. First, the mathematical background of these two estimators and its inference on their respective behavior for long data sets are presented. Second, a simulation study is used to quantify their performance when applied to a reduced number of Doppler data corrupted by noise. A survey of possible improvements of the general mean frequency estimator is then provided, from which an adaptive mean frequency estimator is derived. Finally, application of this estimator to actual Doppler recordings is given.
Since the early stages of ultrasound velocimetry, numerous time-domain estimators of the mean frequency of the Doppler signal have been developed to quantify the mean velocity of the flow. The first estimators were implemented on analogue systems using zero-crossing counters (Rice 1963), (Angelsen and Kristoffersen 1983b) or phase locked loop techniques (Sainz et al. 1976), (Bramanti and Marchesini 1983 ). Thereafter, other signal phase detectors such as the infinite gate detector (Nowicki and Reid 1982), the I/Q algorithm (Barber et al. 1985 ), the double-correlation centroid detector (Reid 1966 ), the single-correlation centroid detector (Arts and Roevros 1972) and, finally, the ~ f centroid detector (Gerzberg and Meindl 1980b) have been used. Recently, estimators with better, but computationally intensive, performances have been developed using autoregressive signal modelling (Kuc and Li 1985); (Loupas and McDicken 1990) or autoregressive moving average modelling of the Doppler signal (Talhami and Kitney 1988). Finally, estimators derived from the properties of the auto-correlation function of the Doppler signal have been developed and applied to color Doppler imaging. Their main advantage consists in offering the best compromise between the quality of mean frequency estimation and the shortness of available signal windows. It must be pointed out that most of the previous estimators can be derived from the proper901
902
Ultrasound in Medicineand Biology
METHODS
and, in particular, for r = 0:
Mathematical background Both estimators CE~ and CE, utilize the properties of the correlation function of a complex value signal. However, the derivation of each mean frequency estimator is based on different assumptions concerning the nature of the signal. Color Doppler estimator CEt The Doppler signal is assumed to be a complex value and continuous time signal: z(t) = x(t) + iy(t)
tE R
( 1)
where R represents the ensemble of real numbers, x(t) and y(t) are the in-phase and quadrature components (I and Q) of the analytic signal, resulting from the baseband conversion of the backscattered radiofrequency signal. This signal can either be assumed deterministic with finite energy, or random with second-order stationarity; subsequently, this difference will have no importance. In all cases, its energy spectrum (deterministic context) or its power spectral density (PSD) (stochastic context), gz(f), is a positive, real value function. The mean frequency of the Doppler signal is defined by:
f=
f
-oD
fg.(f)df
provided that its energy, in a deterministic context, or its mean instantaneous power, in a stochastic context, are normalized:
L
oo
g~(f)df = 1
(3)
~ fg~(f)df
(6)
Therefore, according to eqn (2):
l [d<(,) 1 tTJ,_o
/= E
(7)
Miller and Rochwarger (1972) and Kasai et al. (1985) utilize a modulus and phase representation of
Cz(~): (8)
G ( r ) = I G ( r ) l e '~*) Equation (7) can then be rewritten as follows:
l [ dlG(r)lei't" 1 l[[dlCz(r)l i d~(r) ] ei~(,)] - 2i~ [ -~r + ilG(r) dr J j,=0 (9) Because the magnitude of the Doppler signal IG(OI is an even function, its derivative d] Cz(r) I vanishes dr at r = 0, then:
f:27l[
1
IG(,-)I
l [ ,oU)l
d~ ],=o=2--~[
dr
],=o (lO)
Because [ cf. eqn (3) ]: IC(r)l,=o =
f
--oo
g~(f)df = 1
(11)
o()
of)
Here, E represents the expectation operator. The mean frequency can be determined from the autocorrelation function of the signal, which is the inverse Fourier transform of the PSD:
C=(r) =
= 2i~r
TJ,=o
(2)
c¢3
E{ Iz(t)l 2} =
Volume17, Number9, 1991
g,(f)e2i*Y'df
(4)
Because the phase is an odd function, ¢(r) vanishes at r = 0 so that:
¢(~) lim~o - -
(12)
1 Im{ Cz(z) } lim,-otg- Re{Cz(r) }
(13)
f =~
1
which finally leads to:
c~
Indeed, if the first time-derivative of Cz(r) does exist, we have:
dCz(r) f f o~ d ~ - 2i7r f g z ( f ) eZ"~df oo
(5)
f= ~
l
In this expression, tg -~ represents the arc tangent function. At this stage of the demonstration, one can observe that the C E 1 estimator is perfectly defined for
Adaptive estimation for Doppler signal • A. HERMENTel aL all f G [ - oo, + oo ]. There is no limitation along the range of analyzable frequencies. But, practically, the mean frequency computation must be achieved from discrete samples of the Doppler signal. One can only estimate the correlation function at discrete time instants n = 0, _+1, +_2 • • .. The mean frequency is therefore approximated
903
+oo ~+1/2 ~ C~(n) fe-Zi'f"df n=-oo ~' -1/2
f=
The mean frequency is thus expressed as the sum of the different correlation lags C~(n) weighted by the following coefficients:
by:
+1/2
f f = ~(1)
1
Im{ C~(I)}
(19)
fe-2"f"df=
-1/2
(14)
--
2---~ - 27r tg-I R e { C z ( l ) }
0 for n = 0 (-1).+1 (20) 2i~rn f o r n 4 =0
_
_
which gives: This new expression of f will induce a number of constraints on CEI, which are directly linked to signal sampling. In particular, the range of analyzable frequencies will be reduced to: [ - P R F / 2 , + P R F / 2 ] , where PRF represents the pulse repetition frequency of the Doppler system. Indeed, this first order approximation of the derivative would only be valid provided that the phase is relatively smooth or that it does not change faster than the PRF. Hence, a striking fact about CEI is that the theoretical derivation is based on a continuous signal model, and that eqn (14) is the consequence of the discrete nature of a pulsed Doppler system.
Here, the Doppler signal is assumed to be a complex value and a discrete time signal:
nEZ
(15)
where Z represents the ensemble of integer numbers. The Doppler signal is considered as a sampled version of a continuous time and band-limited signal. Whatever the case (deterministic or random signal), its PSD is a positive, periodic, real value function that can therefore be expanded in a Fourier series:
gz(f) = ~ Cz(n)e -2i'fn
(16)
nEZ
whose Fourier coefficients are defined as:
+l/2
C~(n) =
f -1/2 g=(f)e2~'/"df
nE Z
__
2i7r
X [ ~ -(-l)n+1 Cz(n) + ~ tl~-
oo
.=1
n
(-1)"+1 ] n C=(n) (21)
Taking into account the hermitian symmetry property of the correlation lags of a complex value process:
C=(n) = C*(-n)
(22)
we have:
General mean frequency estimator of CE,
z(n)=x(n)+iy(n)
1 f=
(17)
These coefficients are the correlation lags of the discrete correlation function of the signal. Because g=(f) is periodic, the mean frequency is defined differently, and eqn (2) must be rewritten as follows:
1 ~ ( - 1 ) "+t - f = 25g n
[ C z ( n ) - Cz*(n)]
(23)
and finally: f=l
~ (-l)n+______ ~ Im{Cz(n)}
"/I"n=1
n
(24)
This is an expression of the estimator, previously given by Angelsen and Kristoffersen (1983a). This expression can also be derived from the derivative of the correlation function expressed as a Taylor series expansion. But in spite of this common starting point, the different assumptions made on the nature of the Doppler signal and the way in which the derivative at the origin of C=(T) is estimated in the continuous case, give rise to very different properties for CEI and CE,.
Differences between CEI and CE. Ifa phase and modulus representation is used as in eqn (8), and if N represents the number of available correlation lags, (24) will give:
+zo Cz(n)e-2i'f"df f = ffl5~ fgz(f)df = ~'f+1/2--1/2f,,=~-oo (18)
f=l
~ (-l)"+_____~im{iCz(n)[e,~.)} (25)
71"n=l
n
904
Ultrasound in Medicine and Biology
or;
Volume 17, Number 9, 1991 Cz(n)
1 --
y=
l -
~
(-l)
71" n = l
"+'
IC~(n)l .sin ~o(n)
(26)
1
f = ~ ~o(1)
(27)
Equation (26) will reduce to eqn (27) when: N = 1, ICz(l)l =½
andsin¢(l)=~o(l)
(28)
which is very restrictive. CE~ and CE, are quite different estimators, even if the basic principle of their derivation is the same.
Influence of the correlation lag estimator At this stage the estimators are not completely defined. In practice, the correlation lags are unknown and must be estimated from the Doppler signal samples. The way these correlation lags are computed will influence the properties of the estimator. There are numerous approaches to calculate the correlation lags. However, because of computation time requirements, only two correlation function estimators can be used: 1 N-n
C,(n) ~ E Z,+mZm m=l
forn=0, I...N-
1 (29)
+PRF/2
* Zn+mZm
1 (30)
The expression z* represents the complex conjugate of z. C1 is the conventional correlation estimator that is biased but whose variance is acceptable. The C2 estimator is unbiased. However, its variance becomes high when n is large and close to N. The ratio Im{ C( 1) }/ Re {(C( 1 ) } does not depend on the type of correlation estimator; C~ or C2 can be equally used for CE~. On the other hand, selecting either Ci or C2 for CEn must still be defined. C~ will reduce the contribution oflm{ Cz(n) } with n, thus allowing more rapid convergence of the sum ofeqn (24) and a reduction of the Gibbs phenomenon associated with the truncation. However, this decreasing contribution of Im { (C~(n)} will lead to a bias of the frequency characteristic of CE,, increasing toward the high frequencies. This phenomenon is represented in Fig. 1, which displays the behavior of CE,, for a sine signal with a frequency varying from - P R F / 2 to +PRF/2 and a N = 48 sample window, when the correlation estimators Cl and C2, given in eqns (29) and (30), respectively, are applied. Because CE~ is not dependent on the exact nature of the correlation function estimator and because of the oscillations of the frequency characteristic of CEn associated with the use of the correlation estimator C2, we will subsequently adopt, in this article, the correlation estimator Ci for both the mean frequency estimators. By making this choice, no attention has been paid to the statistical properties (bias and vari-
+PRF/2
J
-PRF/2
Z
forn=0,1..-N-
n
while eqn (14) gives the estimator CEv
N-n
N-n,n= 1
A
i
i
I
+PRF/2
i
-PRF/2
B
•
•
+PRF/2
Fig. 1. Influence of the correlation estimator on CE,,, sine signal, N = 48: x = actual Doppler frequency, y = estimated Doppler frequency. (A) Cj: smooth frequency characteristic, (B) C2 oscillating frequency characteristic.
Adaptive estimation for Doppler signal • A. HERMENTel al.
ance) of the correlation lag estimator. However, the low variance of C~ will be in accordance with the statistical and noisy nature of the Doppler signal.
Theoretical behavior of the estimators An extensive analysis of CEI and CEn has been made by Kristoffersen (1988). The derivation of analytic expressions for the bias and the variance of these estimators in a hypothesis of strong filtering led to the following: CE1 yields a low bias estimate of the mean frequency, which becomes unbiased for symmetrical Doppler spectra. However, the variance of estimation increases rapidly with Doppler signal bandwidth and noise. Finally, these quantities are shown to be independent of the signal mean frequency as far as no aliasing occurs. CE~ is characterized by a bias increasing with the spectrum skewness, the noise and the Doppler signal bandwidth. This bias can reach up to 50% for very large Doppler spectra. However, the noise immunity of the estimator is far better than the one of CEi and, thus, even more when the signal bandwidth is large. Finally, the estimation variance is mean frequency independent as far no aliasing occurs. In addition, Kristoffersen and Angelsen (1985) propose an alternate approach that consists in shifting the spectrum down to the zero frequency by an adaptive premixing of the Doppler signal. This provides an unbiased mean frequency estimator with a variance that is independent of the mean frequency. However, the premixing operation can be applied in a safe way only if the strong filtering hypothesis is valid. This will not always be the case for noisy, wide-band and short data sets of Doppler signal. This is why this adaptive version of CE~ will not be studied in the following developments. However, its performances will be compared to those of other mean frequency estimators on actual Doppler signal recordings. These theoretical results (Kristoffersen 1988) are now illustrated on a synthetic Doppler signal for further comparison with short data sets.
Behavior of the estimators on long data sets Two characteristics are of major importance for a good estimation of the mean frequency of the Doppler signal in color flow mapping: • The frequency characteristic of the estimator, which indicates the range of frequencies for which the relation between the estimated and the true mean frequency of the Doppler signal are reasonably equal, and which is directly linked to the range of flow velocities that the device can manage. • The robustness of the mean frequency estimation
905
in the presence of noise, which directly determines the quality of Doppler images. These parameters are studied for an ideal noiseless signal and for a signal-to-noise ratio ( S / N ) of 0 dB. These values have been chosen from the results of a more extensive study. They are representative of the behavior of estimators when analyzing excellent and rather poor Doppler signals.
Doppler signal model The Doppler signal has been modeled by a complex value signal whose real and imaginary parts are in quadrature to comply with the requirements of the baseband conversion. The Doppler spectrum is rectangular, and its bandwidth can be continuously set between AB = 0 and AB = PRF. This deterministic model of the Doppler signal has been preferred to a stochastic model. It will allow to separate clearly, in the subsequent studies, the respective contribution of the Doppler signal and that of the noise. However, a complementary study, not reported here, and using the Doppler signal model proposed by Mo and Cobbold (1986), has shown that the conclusions provided in this article are comparable for both models. The nature of the noise physically existing on the Doppler signal is hard to model. The real and imaginary components of the electronic noise are mutually uncorrelated. The baseband conversion of the radiofrequency signal leads to a phase difference of 7r/2 between the two components of the noise. Finally, cross-talk phenomena contribute to correlation between these two components. The behavior of the estimators as a function of the noise has thus been studied for white Gaussian and uniformly distributed noises. For each type of noise, equal, in quadrature or completely uncorrelated, real and imaginary parts have been used, and the results obtained in all cases remain identical. The present study corresponds to a white, evenly distributed noise with uncorrelated real and imaginary parts. The frequency characteristics of estimators (mean frequency estimated as a function of the real mean frequency of the Doppler signal) have been obtained by shifting the mean frequency of this signal between - P R F / 2 and + P R F / 2 . The comparison between the behavior of both estimators is achieved for an N = 48 Doppler sample window to eliminate the influence of the signal windowing on the estimation of correlation lags and the truncation effect in the sum of expression (24).
Frequency characteristic Fig. 2 displays the frequency characteristic of
CE~ and that of CEn for a Doppler signal bandwidth
906
Ultrasound in Medicine and Biology
+PRF/2
~'
-PRF/2
Volume 17, Number 9, 1991
+PRF/2
a
i
I
I
A
..o.°o'°
I
+PRF/2
+PRF/2
-PRF/2
•
i
i
i
B
+PRF/2
D
+PRF/2
+PRF/2
°,~oqO°~..~*o°°°
~o ,.°°°°° -PRF/2
C
+PRF/2
-PRF/2
Fig. 2. Frequency characteristics of C E I and CEn:x = actual Doppler frequency, y = estimated Doppler frequency. CEI:dotted line, CE,: plain line. (A) N = 48, AB = 0.001 PRF; (B) N = 48, AB = 0.2 PRF; (C) N = 8, AB = 0.001 PRF; (D) N = 8, AB = 0.2 PRF. Note the inflexions of the frequency characteristic corresponding to _+(PRF-AB)/2, here _+0.4 PRF for CEn in Fig. 2(B).
AB = 0.001 PRF (very narrow band signal) and AB = 0.2 PRF. It clearly shows that CE~ has a frequency range covering the whole Nyquist domain and that CE, has a frequency range that directly depends on the bandwidth of the Doppler signal. These results are in accordance with the theoretical behavior of both estimators (Kristoffersen 1988): The assumption of a continuous time signal does not a priori impose constraints on the frequency range of CEI, which will be able to estimate mean frequencies all along the Nyquist band from - P R F / 2 to + P R F / 2 as far as the phase ~ 1) does not vary too rapidly. In the present case, this estimate is unbiased because the rectangular spectrum chosen is symmetrical with respect to its mean frequency. The assumption of a discrete time Doppler signal and the use of a development in Fourier series of a periodic PSD imply that the mean frequency estimated by CEn will only be significant when the lowest
frequency of the Doppler spectrum will be higher than - P R F / 2 and when its highest frequency will be lower than +PRF/2. Therefore, the frequency range of CE, will be equal to [PRF - AB]. Outside these limits, the bias increases very rapidly along with frequency. Figure 3 presents the average of 64 of the previous frequency characteristics for S / N = 0 dB and AB = 0.2 PRF. One can observe a reduction in the frequency range for both estimators. This can be explained by the increase in the estimation variance when the extreme frequencies of the spectrum cross _+PRF/2. As a result, the averaging of mean frequency, which is commonly used in color flow mapping to decrease the noise on the image, is not really adequate for noisy signals. Indeed, it induces a decrease in the frequency range, which directly depends on the spectral width of the Doppler signal and on the signal-to-noise ratio.
Adaptive estimation for Doppler signal • A. HERMENTel al.
+PRF/2
.
907
+PRF/2
.
.
,
.
-PRF/2
,
•
A
.
+PRF/2
+PRF/2
:
:
;
-PRF/2
.
.
.
.
B
+PRF/2
+PRF/2
.
-PRF/2
:
C
.
.
.
I
+PRF/2
-PRF/2
I
I
.
I
D
.
.
.
+PRF/2
Fig. 3. Average m e a n frequency characteristic of CE, and CEn for a noisy Doppler signal, S / N = 0 dB, AB = 0.2 PRF: x = actual Doppler frequency, y = averaged estimated Doppler frequency. ( A ) CE,, N = 48; (B) CE~, N = 48; ( C ) CE~, N = 8; ( D ) CE,, N = 8.
Estimation variance Figure 4 presents the standard deviation (SD) of a mean frequency estimated as a function of the signal mean frequency, for AB = 0.2 PRF and S / N = 0 dB. Figure 5 presents the SD for these two estimators as a function of the signal bandwidth for S / N = 0 dB, but for a fixed zero mean frequency. Again, these results are in accordance with the theoretical behavior of both estimators (Kristoffersen 1988), in particular, the rather constant variance of CE1 as far as no aliasing occurs, and the strong increase in this variance when the extreme frequencies of the Doppler spectrum cross the Nyquist limits. CEn has a far better noise immunity than CEI, and this is even more observed as the signal bandwidth increases. The performances of the estimators directly depend on their ability to estimate the derivative at the origin of the correlation function. If the correlation function varies slowly around the origin, the correla-
tion lags evolve slowly. The derivative is correctly estimated from C( 1 ) alone for CEI and from the sum of eqn (24) for CErt. If the variations of the correlation function at the origin become faster (high mean frequency or wideband signal), the derivative will be less and less well estimated by CEI. On the other hand, the correlation lags decrease faster, the sum ofeqn (24) will converge faster, and CEn will adequately estimate the mean frequency in the [ - ( P R F - AB)/2, + ( P R F - AB)/2] range, thus accounting for a lower variance. In addition, eqn (14) shows that CE~ only uses the first correlation lag, and any uncertainty on the value of Cz( 1 ) is directly converted into an uncertainty on the mean frequency. The impact of this uncertainty will be minimal when Im{ Cz( 1 )} will be small and when Re{ Cz(1)} will be large, i.e., to say for the low mean frequencies; conversely, it will be maximal around __+PRF/2. On the other hand, CEn utilizes a weighted aver-
908
Ultrasound in Medicine and Biology +PRF/2
Volume17, Number 9, 1991 +PRF/2
'
lr
i
-PRF/2
A
+PRF/2
+PRF/2
:
'
I
e
I
I
B
I
+PRF/2
+PRF/2
., -PRF/2
'
;
-PRF/2
C
+PRF/2
-PRF/2
:
_,
:
:
D
;
:
:
+PRF/2
Fig. 4. Standard deviation of CEI and CE, as a function of the signal mean frequency, S/N = 0 dB, AB = 0.2: x = actual Doppler frequency, y = standard deviation of the estimated Doppler frequency. (A) CE~, N = 48; (B) CE,, N = 48; (C) CE~, N = 8; (D) CE,, N = 8.
age of the Im{ Cz(n) }, and the uncertainties in calculation o f the Cz (n) may be expected to decrease in the summation of eqn (24).
Influence of the truncation of the Doppler signal T o comply with real-time color Doppler imaging requirements, the n u m b e r of samples available to calculate the mean frequency is reduced ( 4 - 1 2 samples). In the following, if not otherwise mentioned, Nwill be taken equal to 8 for short windows. This will have a direct impact on the sum in eqn (24) for CE, and on the estimation of correlation lags for both estimators.
Influence of truncation on the mean frequency expression The truncation on the sum in eqn (24) will lead to an increased Gibbs phenomenon, which will generate oscillations in the frequency characteristic, even if the correlation lags are perfectly estimated. This phenomenon is illustrated in Fig. 6. It displays the fre-
quency characteristic of CE, for the narrow band signal (zkB = 0.001 ) and when an N = 256 window is used for correlation lags computation. Figures 6 (A,B) corresponds to the utilization o f the only first 48 and the only first 8 of these previous correlation lags in eqn (24), respectively. The oscillations in the frequency characteristic are even more important as the number of lags is low.
Influence of the window length on the frequency characteristic Figure 2 allows one to compare the frequency characteristics of CE, and that of CEn for long and short windows for 2xB -- 0.001 PRF and a d / = 0.2 PRF. Comparison of these curves shows that the frequency characteristic of CE1 is independent of the window length and that the one o f CE, depends on the window length, besides the already observed dependence of this frequency characteristic with the signal bandwidth for CE,.
Adaptive estimation for Doppler signal • A. HERMENTet aL
909
0,3
// 0,2
0,1
B'----~
0,0
0,0
r-
0,2
0,4
0,6
0,8
1,0
Fig. 5. Standard deviation of the mean frequency estimation by CE~ and CE, for a noisy signal as a function of spectrum width; S/N = 0 dB: x = actual Doppler bandwidth, y = standard deviation of the estimated Doppler frequency. (A) CE,, N = 48; (B) CE,, N = 48; (C) CE~, N = 8; (D) CE,, N = 8.
Figure 3 compares the average m e a n frequency characteristics of CE~ and CE, for long and short windows, for AB = 0.2 P R F and for a 0 dB signal-to-noise ratio. The general influence of the decrease in the window length is a m o r e i m p o r t a n t degradation of the average frequency characteristic for CE~ than for CE,, even though the frequency range o f CE~ is always better than the one o f CE, for a given window.
+PRF/2
Influence of the window length on the estimation variance Figure 4 indicates the modification o f the SD o f m e a n frequency estimation o f CE~ and CE, along with the frequency, for long and short windows, a n d for S / N = 0 dB. For both estimators, the SD increases with the reduction of the window length. For longer windows, the correlation function o f
+PRF/2
f
/
-PRF/2
A
+PRF/2
-PRF/2
B
+PRF/2
Fig. 6. Influence of a truncation in eq. (24), N = 256, AB = 0.001 PRF: x = actual Doppler frequency, y = estimated Doppler frequency. (A) 48 correlation lags; (B) 8 correlation lags.
910
Ultrasound in Medicine and Biology
the noisy part of the signal tends toward the theoretical correlation function of the noise, which has no contribution on the correlation lags excepted on Cz(0). In addition, the SD of CEn remains rather constant with frequency, whatever the window when that of CE~ is degraded around + / - PRF/2 for short windows as a result of the CE~ behavior itself and of an artificial enlargement of AB caused by window truncation. Figure 5 compares the SD of these estimators for the same 0 dB signal-to-noise ratio and the same N = 8 window. Again, one can observe that the SD increases with the decrease of N, CE~ being more sensitive to the window reduction. Moreover, the SD remains independent of the signal bandwidth for CEn, while it is extremely AB-dependent, tending toward a constant value when the signal bandwidth increases, for CE1, whatever the N value. Starting from the results of this comparison, the question of defining an estimator allowing to conciliate a high-frequency range and a low-estimation variance is now raised. For that purpose, the starting point will be CEn given by eqn (24), CE~ being completely defined by eqn (13).
Optimization of CEn Three solutions have been studied to improve the
CEn estimator: first, to replace the weighting coefficients an = (-l)"+~/a'n by coefficients taking into account the Gibbs phenomenon associated with the expansion truncation; second, to extrapolate the correlation function so as to create a greater number of correlation lags; and, third, to interpolate the correlation function for the same purpose.
Modification of a, coefficients Here, the problem consists in finding the best polynomial approximation of the ideal frequency characteristic. It has been shown that using the coefficients an = (-1 )n + l / r n associated with the unbiased correlation estimator C2 of eqn (30) led to a larger frequency range but to a more oscillatin~ characteristic of the estimator. In fact, changing C~ to C2 gives the same estimator equation as changing the theoretical weighting sequence { an } n = 1, 2 , . . . , N - 1 for the following:
Volume 17, Number 9, 1991
prove the frequency range and increase the Gibbs phenomenon.
Exponential weighting A simple means to reinforce the frequency characteristic of the estimator is to choose a weighting sequence such as: ( - 1 ) n+,
an -
an
( - 1 ) n+l
N- n
lrn
(32)
with ~ < 1. Figure 7 displays, for N = 8, for AB = 0.001 PRF, and AB = 0.4 PRF, a comparison between the estimator CEn of eqn (24) and the exponentially weighted estimator defined in eqn (32) for/3 = 0.6. The absolute maximum value of the estimated mean frequency increases along with the decrease of /3, but the frequency range remains unchanged. The frequency range of the new estimator remains, in fact, directly connected to the bandpass of the Doppler signal as for CEn.
Least squares weighting A more subtle approach has been proposed by Angelsen and Kristoffersen (1983a). It consists in determining the series of coefficients dn so as to minimize the mean quadratic error on:
J,=o that is to say, to calculate the dn coefficients minimizing:
LI L d.,- J,.=o L d,-J,-oolJ
(33)
where E represents the expectation operator, Cz(z) the estimated correlation function and Cz(r) the true correlation function. For that purpose, one must define a probability law for the PSD obtained in order to calculate the previous mean quadratic error:
---&---, j,=o - t---&-, J,=o =
1
71-na
Dk(f)gz(f)df
(34)
(31)
while keeping the same C~ estimator. One can argue that any change in the sequence { an } that will increase the relative weight values of an with n, will i m -
where G(f) represents the frequency band of the Doppler signal and where: k
Dk(f) = X [d,e -z'f'] - 2iTrf n=-k
(35)
Adaptive estimation for Doppler signal • A. HERMENTet al.
+PRF/2
.
911
+PRF/2
.
.
i
.
-PRF/2
;
A
;
I
+PRF/2
+PRF/2
i
i
I
i
-PRF/2
I
I
B
I
+PRFI2
+PRF/2
I
-PRFI2
C
+PRF/2
-PRF/2
I
I
D
I
+PRF/2
Fig. 7. Influence of an exponential weighting of the a, series on CEn, N = 8. x = actual Doppler frequency, y = estimated Doppler frequency. (A) A B = 0.001 PRF;/~ = I;(B) A B = 0.001 PRF, ~ = 0.6;(C) A B = 0.4 PRF;/~ = l; (D) A B = 0.4 PRF, B = 0.6.
By choosing an even probability law on G(f) and a line spectrum gz(f) = ~( f ) , the minimizing criterion is then written as follows:
min
d,, = fc
IDk(f)i2df
+PRF/2
A
,J"
(36)
•/
\
(f)
Therefore, the dn are the coefficients o f the trigonometric polynomial that is closest, to 2 h r f o n G ( f ) , in the least squares sense. One can note that if G(f) = [ - P R F / 2 , + P R F / 2 ] , one again finds the coefficients an = (-1)"+1/2mr, which are the Euler coefficients for expanding 2i~-fin a Taylor series. In this approach, if G(f) reduces to [ - f 0 , + f 0 ] , the approximation will be better on G(f) iff0 is low for a given n. Figure 8 displays a comparison, for AB = 0.001 PRF, N = 8 samples and for G(f) = [ - 0 . 4 5 PRF, +0.45 P R F ] , between the frequency characteristics o f CE, in dotted lines and the one corresponding to the
j,"
f
j
"~
t"' '/ ."
-PRF/2
~
I
|
0
I
!
I
+PRF/2
Fig. 8. Influence of a least squares optimization of the an series, N = 8, AB = 0.001 PRF: x = actual Doppler frequency, y = estimated Doppler frequency. (A) an = ( - 1) n+, / ~rn; (B) an coefficients according to Angelsen and Kristoffersen (1983a).
912
Ultrasound in Medicine and Biology
coefficients given by Angelsen and Kristoffersen (1983a), as a broken line. The optimized characteristic of the estimator is more linear, but its frequency range is still unsatisfactory. This approach does not either seem to Satisfy the requirements as originally stated. Kristoffersen (1988) applied the Same least squares approach to design a mean frequency estimator offering a constant ~f/ferror characteristics. The result of this approach is comparable to the previous exponential weighting: the maximum value of the curve increases, but its frequency abscissa remains the same.
Extrapolation of the correlation function for n > N Signal modeling allows one to extrapolate the correlation function for n > N. The most direct method (Kay and Marple 1981 ) models the correlation function Cz(n) by an AR model of adequate order, and extrapolates Cz(n) by a Burg algorithm that allows one to define, among all the possible extrapolations of the correlation function, the one that has a maximum entropy (Schwindlein and Evans 1989). Beyond N, however, extrapolation must be limited because the stability of estimation cannot always be ensured. A good compromise would consist of taking, for Doppler signal analysis, [0 < n < 2N] (Schwindlein and Evans 1989). However, the complexity of this approach does not bring a significant gain on the estimation of mean frequency. The extrapolated values of the correlation
function for n > N have a weak ( 1/~rn) contribution in eqn ( 24 ). Last, the order (P) of AR models must be limited to P < N for an estimation of this model by least squares techniques. For Doppler signals, it is generally impossible to obtain a ratio higher than P / N = 1/4 (Kaluzinski 1989). One must notice that color flow mapping leads to 4 > N > 12, i.e., AR model orders smaller than 4, which do not allow for a correct Doppler signal model. This model order limitation can be removed with long AR models (P = N) associated with regularized least squares methods (Demoment et al. 1988). The regularization principle, applied here, consisted in estimating the AR parameters while minimizing a distance that not only depends on the fidelity of the model to the data (least squares approach), but also on a quadratic criterion characterizing the AR spectrum smoothness. The extrapolated series of correlation lags is then obtained from AR coefficients using an inverse Levinson algorithm. Despite this additional effort, the results obtained are still almost equivalent to those directly provided by eqn (24). The frequency range is still limited by the window length as shown in Fig. 9 for AB = 0.001 PRF and AB = 0.4 PRF. In the present case, it is likely that the increase in the model order contributes in itself to a better description of the Doppler signal but that, on the other hand, the AR coefficients obtained with regularization techniques are mutually correlated and that the overall result remains equivalent to that of the previous approach. To conclude, the complexity introduced by these methods does not seem to be justified with respect to +PRF/2
+PRF/2
-PRF/2
Volume 17, Number 9, 1991
A
+PRF/2
-PRF/2
B
+PRF/2
Fig. 9. Correlation lag estimation by AR modelling, N = 8: x = actual Doppler frequency, y = estimated Doppler frequency. (A) AB = 0.001 PRF; (B) AB = 0.4 PRF.
Adaptive estimation for Doppler signal Q A. HERMENT el al.
913
the possible improvement it would bring to the mean frequency estimation.
ues C~(n) = C~(3m) by writing Cz(n) again in the following form:
Interpolation of correlation lags An interpolation of the correlation lags C~(n) will not only lead to more available correlation lags to estimate f, but also to create new intermediate correlation lags corresponding to a higher sampling frequency, because their delay is a submultiple of the sampling period. This will not eliminate aliasing problems beyond PRF/2, but it will enable the operation on the linear part of the frequency characteristic. For instance, creating a correlation lag interpolated between each correlation lag Cz(n) would lead to an equivalent N = 16 window and only use the f values from - P R F / 4 to + P R F / 4 on the frequency characteristic. However, as the vertical scale is doubled because of oversampling, this will correspond to a real range of - P R F / 2 to + P R F / 2 while keeping the linearity of the center part of the frequency characteristic of CE, between - P R F / 4 and + P R F / 4 . The first option consists in defining the number of correlation lags to interpolate between successive original correlation lags C, (n). Trials have been made on different orders. It appears that the addition of two lags interpolated between two consecutive Cz(n) allows, for short windows, optimization of the ratio between the quality of estimation of the mean frequency and the amount of additional calculations introduced by interpolation. The second option consists in determining the interpolation scheme. To obtain an estimator with a high-frequency range, we have chosen to elaborate a scheme allowing to get closer to hypotheses underlying to the setting up of estimator CE~, i.e., an interpolation that tends to reduce the variation of the modulus of the correlation function at the origin and that strengthens the influence to the phase of the first correlation lag. For that purpose, the modulus and the phase of the interpolated correlation lag were linearly interpolated between the successive original correlation lags. One would note that this interpolation is fundamentally different from a linear interpolation of the real and imaginary parts of the Cz(n). The interpolated estimator CE3n is thus defined by:
Cz(n) = ICz(n)le i~")
2
f = l -- N-| E [(-l)3m-t+la3m-tlm{C'~(3 m - l ) } ] 71" m=l /=0
(37) The values ofC'( 3m - l) are calculated from the val-
(38)
defined by its modulus ICz(n) I and its argument ~o(n). The interpolation scheme is: IC'z(3m)l = ICz(n)l
(39)
I C ~ ( 3 m - 1)1 = [2[Cz(n)l + I C z ( n - 1)1]/3
(40)
=[[Cz(n)l+2lCz(n-1)l]/3
(41)
C'z(3m - 2)1
~0'(3m) = ~o(n)
(42)
~'(3m - 1 ) = [2~o(n) + ~o(n - 1)1/3
(43)
~o'(3m - 2) = [~0(n) + 2~o(n - 1)]/3
(44)
C'z(3m) = IC'~(3m)le i¢~3m)
(45)
C'z(3m- 1 ) = IC'z(3m- 1)l e~¢t3'~-l)
(46)
C'z(3m - 2) = IC~(3m - 2)le i¢~3"-2)
(47)
The a3m_l coefficients are also interpolated by: aam = 3 / 3 m 1)
(49)
a3m_2 = 3 / ( 3 m - 2)
(50)
a3m_ 1 =
3/(3m-
(48)
The frequency characteristic of CE3n is identical to that of CE~ for AB < 0.4 and slightly differs for wider spectra. However, as illustrated in Fig. 10, even for large bandwidth signals (AB = 0.9 PRF) and a short window, there is an obvious similarity between the frequency characteristics of both estimators. As for the estimation variance, it is strictly equivalent for both CE3n and CE~. Finally, note that ifa3m_l = a3m-2 = 0 , then CE3n = CE~. Hybrid estimator Equations (37-50) allow for a transition from estimator CEn to estimator CE1. Therefore, it seems possible to define, by an appropriate weighting of coefficients a3m-1 of eqns (48-50), an estimator halfway between CEn and CE~ and, in particular, to adapt the compromise between the frequency range and the variance of estimation to the characteristics of the sig-
914
Ultrasound in Medicine and Biology
+PRF/2
+PRF/2
-PRF/2
Fig. 10. Comparison between CE, and C&, N = 8, AB = 0.9 PRF: x = actual Doppler frequency, y = estimated Doppler frequency. Dotted line: CE, ; plain line: CE,,
nal to process. We have thus generated a hybrid estimator CE,,, defined as follows:
f=
b2;[(-1)3”-‘+‘aJ,_,(cX)
Volume 17, Number 9, 199 I
However, if some characteristics on the Doppler equipment, such as the noise level of the system are known a priori, or if some information on ‘the Doppler signal characteristics can be inferred from the environment of the considered window, the characteristics of the mean frequency estimator can be adapted, for each window, to process the corresponding signal. The following scheme presents one simple application of this adaptive estimator to the Doppler signal. The adaptivity is based on the assumption that the Doppler signal mean frequency cannot be very different between two adjacent windows:
.L+,%fk
(55)
In the case of a monogate or TM Doppler system, k will characterize two windows located at the same site but at consecutive instants. For a color Doppler imaging system, k will correspond to two windows on different sites: either two consecutive windows along the same scanning line, or two windows located at the same depth along two consecutive scanning lines. The criterion selected to simply express the gradual changes of the Doppler signal is as follows:
I=0 X Im{C;(3m
-
f)j]
(51)
The values of C:( 3m - I) are calculated with eqns (38,47), but a3m_1coefficients are replaced by new coefficients u3,,_,( LY)defined as follows: ~~~(a) = (1 +
2&)/3m
(52)
a3m-1(cu) = a( 1 + 2a)/( 3m - 1)
(53)
= (Y(1 + 2a)/(3m - 2)
(54)
a3&~)
(Yis a trimming parameter ranging between [ 0 and l] and determined by the user.
Adaptive estimator An improvement can thus be obtained for a given signal by choosing the value (Ythat is most appropriate to the characteristics of this signal. Unfortunately, the statistical properties of the Doppler signal change very rapidly throughout the cardiac cycle, even for a given location of the flow. The signal-tonoise ratio strongly varies between diastole and systole as a result of hemodynamical phenomena inherent in flow conditions and to the echo canceller cutoff characteristics. Moreover, the signal bandwidth can also be strongly modified according to flow conditions during the same heart contraction, which leads to a difficult choice for (Y.
(56) a2 allows one to smooth the variation of fbetween Sk and Jik+,. A low value for a2 impedes any rapid modification of the mean frequency between two consecutive signal windows. As for the coefficient QI, , it avoids an oversized CJ (Y, , a2) in eqn (56) when the mean frequency is close to zero and permits mean frequency variations when fk = 0. An illustration of the respective contribution of these two parameters is given in Fig. 11. The detailed block diagram of Fig. 12 displays an algorithm that chooses the largest (Yvalue complying with the smoothness criterion Cr( a,, a2). a0 represents the incremental step of CL
al
t&small
a2 too small
Fig. 11. Influence of a, and (Yeon the mean frequency estimation. x = time; y = estimated Doppler frequency.
Adaptiveestimationfor Dopplersignal• A. HERMENTet al.
915
+PRF/2
k=0 f0=0 OCO,(gl, 0[2
a=l+ao
k=k+l
I
k-
•
I
l
I
•
•
•
+PRF/2
Fig. 14. Average mean frequency characteristic of CE3,,,,, CE~ and CE,, for a noisy Doppler signal, S / N = 0 dB, AB = 0.2 PRF, N = 8: x = actual Doppler frequency, y = averaged estimated Doppler frequency. (A) CEI; (B)
f = EC3n~
fEn;
F,
I
-PRF/2
(C) CE3na,
ot =
0.8.
Cf ( a l , tI2)
lno
~=0
RESULTS
I no
~ i Y es k=k~
I no
yes
e
END
I
Fig. 12. Principle of adaptive mean frequency estimation with CE3n~.
Properties of the hybrid estimator Figures 13-16 display the b e h a v i o r o f CE3n . for a = 0.8, N = 8 a n d AB = 0.2 P R F together with those o f CE~ a n d CE, for c o m p a r i s o n o f their respective performances. Analysis o f these curves does show, as expected, for CE3n,~, an intermediate b e h a v i o r between that o f CE t a n d that o f CE, b o t h for the bias a n d the variance. A n y a value larger t h a n 0.8 will give to CE3,,, a behavior closer to that o f CEI a n d a n y a value lower t h a n 0.8 a behavior closer to that o f CE,.
+PRF/2 +PRF/2
I
-PRF/2
+PRF/2
Fig. 13. Frequency characteristics of CE3n,,, CEt and CE,,, N = 8, AB = 0.2; PRF: x = actual Doppler frequency, y = estimated Doppler frequency. (A) CE~; (B) CE,,; (C) CE3.,,, a = 0.8.
-PRF/2
I
I
I
+PRF/2
Fig. 15. Standard deviation of C3,~, CEi and CE, as a function of the signal mean frequency, S / N = 0 dB, AB = 0.2 PRF, N = 8: x = actual Doppler frequency, y = standard deviation of the estimated Doppler frequency. (A) CE~ ; (B) CE.; (C) CE3,,,,, a = 0.8.
916
Ultrasound in Medicine and Biology
0,2
Volume 17, Number 9, 1991
0,6
094
098
190
Fig. 16. Standard deviation of the mean frequency estimation of CE,,,. CE, and CE, for a noisy signal as a function of spectrum width; N = 8, S/N = 0 dB: x = actual Doppler bandwidth, y = standard deviation of the estimated Doppler frequency. (A) CE, ; (B) CE,,; (C) CE,,,, (Y= 0.8.
Application of the adaptive estimator to actual Doppler signals Doppler signal recordings have been achieved using a pulsed velocimeter operating at 3 MHz. Two sample volumes were located downstream from a stenosis, one in the laminar jet and the other in the area of turbulent flow corresponding to the breakup zone
of the jet. These recordings have been made in vitro using a pulsatile flow. In the following application, N = 4 sample windows were used to test the adaptive estimator under stringent conditions. The gap between two consecutive windows k and k + 1 was of 60 sampling periods as explained in Fig. 17.
1--11111-111-11
N=60
N=4
Doppler window Tk
N=60
Unprocessed data
N=4
N=60
Doppler wicidow fk+l
Fig. 17. Timing used for adaptive processing of actual Doppler signals. x = time, y = Doppler signal amplitude.
Adaptive estimation for Doppler signal • A. HERMENTel al.
+PRF/2
917
+PRF/2
1 •"
-PRF/2
|
l
A
0.75 s
+PRF/2
-PRF/2
B
+0.75 s
+PRF/2
'I
-PRF/2
C
0.75 s
t-PRF/2
|
0/ .1
ii
|,,,
D
0.75 s
+PRF/2
Fig. 18. Behavior of CE~, CE, and CE3,,, on Doppler signals: laminar flow. x = actual Doppler frequency, y = estimated Doppler frequency. (A) CE~, N = 48; (B) CEt, N = 4; (C) CE,,, N = 4; (D) CEa,,,,, N = 4; (E) CE,,, premixing, N = 4.
-PRF/2
E
0.75 s
The results, obtained using the algorithm described in Fig. 12, are compared with those of C E i and CE,, for laminar flows on Fig. 18 and for a turbulent flow on Fig. 19. On these two figures, curves A correspond to the mean frequency estimation using CE~ and an N = 48 sample window. They will be considered as the exact mean frequency of the Doppler signal even if this is a rough approximation
in the case of a turbulent flow. Curves B, C, D and E correspond to CE~ , CE,, CE3n a and CE,,, when using the premixing scheme proposed by Kristoffersen (1985). In use of the adaptive estimator, the same values at = 0.001 and a2 = 0.5 were used for both laminar and turbulent flows. In this premixing scheme, the frequency shift of the Doppler signal was simply taken
918
Ultrasound in Medicine and Biology +PRF/2
Volume 17, Number 9, 1991 +PRF/2
r-
""]
-PRF/2
A
r
0,75 s
I'!t -PRF/2
I' lpl B
+0.75 s
D
0.75 s
+PRF/2
+PRF/2
-PRF/2
C
0.75 s
-PRF/2
+3 PRF/2
Fig. 19. Behavior of CEI, CE, and CE3,~ on Doppler signals: turbulent flow. x = actual Doppler frequency, y = estimated Doppler frequency. (A) CE~, N = 48; (B) CE,, N = 4; (C) CE,,, N = 4; (D) CE3,,a,N = 4. (E) CE,,, premixing, N = 4.
O...l
~r'-~ -PRF/2
I E
0.75 s
equal to the mean frequency fk. estimated in the previous time window. The adaptivity of the CE3n,~estimator allows one to improve significantly the quality of estimation of the mean frequency by adequately estimating the highest frequencies of the systolic phase while keeping a reduced estimation variance during the diastole.
When the signal premixing down to the zero frequency is applied to the laminar flow Doppler signal, the mean frequency is excellent. However, for the turbulent flow, an important mean frequency error propagates all along the cycle because o f the transmission of local erroneous mean frequency estimations. From these results, improvement of color flow
Adaptive estimation for Doppler signal • A. HERMENTel aL
mapping images can be contemplated. Adaptive mean frequency estimation should improve the image quality for a given frame rate, but also increase the frame rate while keeping comparable images. This should be of some interest for cardiac imaging where the flow modifications are extremely fast during the systole. Finally, it can offer a partial contribution to improvement in the compromise between low-frequency cutoff and rise time in the echo cancellers. The excellent noise immunity of the adaptive estimator for low frequencies, as for the CE,,, should allow better analysis of low frequencies. CONCLUSION The color Doppler mean frequency estimator
CE~ and the general mean frequency estimator CE, are based on the evaluation of the first derivative of the Doppler signal correlation function at the origin. The first one has a frequency range of [ - P R F / 2 , + P R F / 2 ], but its variance of estimation is important and its noise immunity very poor. The second one has a reduced frequency range but a low estimation variance and an attractive noise immunity. Using a nonlinear interpolation of the different correlation lags for the Doppler signal correlation function derivative estimation, has led to a new mean frequency estimator. Its behavior is continuously shifted from that of the CE~to that of the CE,, according to the value of a single trimming parameter. This hybrid estimator allows one to optimize, for a given window, the compromise between the frequency range and the estimation variance to the particular features of a given Doppler signal. Its preliminary adaptive application to the Doppler signals reveals its interest when one is processing very short Doppler windows. However, the best adaptivity criterion has still to be defined for each particular application, in particular, adaptivity schemes fitted to color Doppler imaging. Acknowledgments--We gratefully acknowledge Mrs. F. Marchand for her editorial and secretarial assistance.
REFERENCES Angelsen, B. A.; Kristoffersen, K. Discrete time estimation of the mean Doppler frequency in ultrasonic blood velocity measurements. IEEE Trans. Biomed. Engin. BME-30:207-214; 1983a. Angelsen, B. A.; Kristoffersen, K. Zero-crossing density for ultrasonic Doppler signals obtained from computer simulations. Ultrasound Med Biol 9:661-665; 1983b.
919
Arts, A. G. J.; Roevros, J. M. On the instantaneous measurement of blood flow by ultrasonic means. Med. Biol. Eng. 10:23-34; 1972. Barber, W. D.; Eberhard, J. W.; Carr, S. G. A new time-domain technique for velocity measurements using Doppler ultrasound. IEEE Trans. Biomed. Eng. BME-32:213-229; 1985. Bramanti, M.; Marchesini, E. Analysis ofa PLL system as analog processor in CW Doppler flowmeters. IEEE Trans. Biomed. Eng. BME-3:584-589; 1983. Demoment, G.; Houacine, A.; Herment, A.; Mouttapa, I. Adaptive Bayesian spectrum estimation. Fourth ASSP Workshop on spectrum estimation and modeling, Minneapolis, USA, August 3-5; 1988. Gerzberg, L.; Meindl, J. D. Power spectrum centroid detection for Doppler systems applications. Ultrason. Imaging 2:232-261; 1980a. Gerzberg, L.; Meindl, J. D. The ~ power spectrum centroid detector: System considerations, implementation and performance. Ultrason. Imaging 2:262-289; 1980b. Kaluzinski, K. Order selection in Doppler blood flowsignal spectral analysis using autoregressive modelling. Med. Biol. Eng. Comput. 27:89-92; 1989. Kasai, C.; Namekawa, K.; Kouano, A.; Omoto, R. Real-time twodimensional blood tlow imaging using an autocorrelation technique. IEEE Trans. Sorties Ultrasonics SU-32:458-464; 1985. Kay, S.; Marple, S. L. Spectrum analysis. A modern perspective. Proc. IEEE 69:1380-1419; 1981. Kristoffersen, K. Time domain estimation of the center frequency and spread of the Doppler spectra in diagnostic ultrasound. IEEE Trans. Ultrason. Ferroelec. Freq. Contr 35:685-700; 1988. Kristoffersen, K.; Angelsen, B. A. J. A comparison between mean frequency estimators for multigated Doppler systems with serial signal processing. IEEE Trans. Biomed. Eng. BME 32:645-657; 1985. Kuc, R.; Li, H. Reduced-order autoregressive modeling for centerfrequency estimation. Ultrason. Imaging 7:244-251; 1985. Loupas, T.; McDicken, N. Low order complex AR models for mean and maximum frequency estimation in the context of Doppler color flow mapping. IEEE Trans. Ultrason. Ferroelec. Freq. Contr 37:590--601; 1980. Miller, K. S.; Rochwarger, M. M. A covariance approach to spectral moment estimation. IEEE Trans. Information Theory 18:588596; 1972. Mo, L. Y.; Cobbold, R. S. Speckle in continuous wave Doppler ultrasound spectra: A simulation study. IEEE Trans. Ultrason. Ferroelec. Freq. Cont. 33:747-753; 1986. Nowicki, A.; Reid, J. M. An infinite gate pulsed Doppler. Ultrasound Med. Biol. 7:41-50; 1981. Reid, J. M. Processing and display techniques for Doppler flow signals. In: Harrison, D. C.; Sandier, H.; Miller, H. A., eds. Cardiovascular image and image processing. The society of photooptical instrumentationEngineers; 1966:73-78. Rice, S. O. Mathematical analysis of random noise. BSTJ 23:288; 1963. Sainz, A.; Roberts, V. C.; Pinardi, G. Phase-locked loop techniques applied to ultrasound Doppler signal processing. Ultrasonics 14:128-132; 1976. Schlindwein, F. S.; Evans, D. H. A real-time autoregressive spectrum analysis for Doppler ultrasound signals. Ultrasound Med. Biol. 15:263-272; 1989. Talhami, H. E.; Kitney, R. I. Maximum likelihood frequency tracking of the audio pulsed Doppler ultrasound using a Kalman filter. Ultrasound Med. Biol. 14:599-609; 1988.