System for Estimating Target Velocity From Pulse Echoes in Response to Their Correspondence With Predetermined Delay Trajectories Corresponding to Different Distinct Velocities

The present invention relates to systems (methods and apparatus) for estimating the velocity of moving targets and particularly the velocity of internal body structures such as tissue and particularly moving blood (blood flow). The invention utilizes the principles of pulse echo systems wherein bursts of signals are repetitively transmitted, propagate to a target and echoes or returns from the target are received, and particularly in which time domain processing of the returns is utilized, rather than frequency domain processing in which the frequency shift due to the Doppler effect is obtained.

In accordance with the invention, however, the returns are processed to search return echo signals for predetermined target trajectories. Conceptually in two dimensional space defined by time coordinates, one being the discrete, ordered times of transmission and the other the times of return of the echoes or the range of the target, the sets of times are along different realizations of a plurality of lines of different slope, each corresponding to a different target velocity. The lines are, thus, along different delay trajectories. The line or delay trajectory along which a plurality of successive returns is most proximal is found by searching over a plurality of different delay trajectories of the returns to determine the most likely correspondence to the expected return from a target traveling along that trajectory. The velocity corresponding to that trajectory provides an estimate of the velocity of the target.

The invention is especially suitable for use in diagnostic ultrasonic imaging and/or blood flow measurement instruments. The invention, however, is generally useful wherever pulse echo signals are transmitted and returns received and wherever the measurement of the velocity of a target is of interest.

Doppler techniques involving Doppler frequency shifts and time domain measurements on echoes have been used for velocity determination in medical diagnostic instruments, particularly where blood flow velocity is of interest. The noisy, clutter full environment in which such techniques are operative makes their application difficult. A system utilizing such techniques which has noise immunity and is robust enough to provide reliable accurate velocity estimates has long been, and is still, sought after by workers in the field. It is a feature of this invention to provide such a system.

Pulsed Doppler systems use the change of phase from successive reflections to estimate the velocity of the moving scatterer. Peak velocity can be found using the Fast Fourier Transform (FFT). Other methods like autoregressive spectral estimation can be used to find the peak velocity. A review of the Doppler frequency shift techniques for velocity estimation may be found in an article by Vaitkis, et al, "A Comparative Study and Assessment of Doppler Ultrasound Spectral Testimation Techniques, Pt. 1, Estimation Methods, in Ultras. Med. & Biol., 14, 661-672 (1988). Mean velocity can be found using autocorrelation techniques disclosed in C. Kasai et al., "Real-time 2-Dimensional Blood Flow Imaging Using an Autocorrelation Technique," IEEE Trans., on Sonics & Ultras., SU-32, 458-463 (1985). A representative implementation of such techniques may be found in Ishikawa, Katakura, and Ogawa, U.S. Pat. No. 4,809,703, issued Mar. 7, 1989. A different set of velocity estimators seek to determine the time shifts between successive transmission, using cross-correlation techniques. Time domain techniques have been reviewed in an article by Hein and O'Brien. " Current Time Domain Methods for Assessing Tissue Motion by Analysis from Reflected Ultrasound Echoes--A Review," IEEE Trans. UFFC, Vol. 40, pages 1-19 (1993). The implementation of such techniques in practical forms which are designed to accommodate noise and provide reasonable resolution, say of the order of a centimeter per second (cm/sec), have needed complex computational processes, such as correlations or transforms, in order to carry out the algorithms on which they are based. The problem therefore remains of providing a robust, noise immune system, which may be implemented with less complex and/or more efficient computational software and hardware than heretofore believed to be necessary by workers in this field, and it is a feature of this invention to provide an improved system for processing pulse Doppler returns which is more robust, and less complex than known systems for the same purpose and can operate in a noisy environment so that the improved system is a practical solution for blood flow velocity measurement.

The pulse Doppler measurements which are made in accordance with the conventional techniques, compute target velocity either from the spectrum of the returns or the time of occurrence of the returns at a specific range. The range is determined by the round trip delay time to the point of transmission of the ultrasonic RF pulses, usually the location of the transducer which both transmits the pulses and receives the return. In the case of frequency measurement techniques, usually using the FFT as a spectrum estimator, the maximum or mean single frequency in the spectrum at the preset range is determined. Turbulent flow is considered to have many velocities and may spread the spectrum. Variances of the measured velocities are computed in accordance with specified algorithms.

In contrast to these known techniques, the invention estimates velocity by finding upon which of a plurality of velocity lines or trajectories, each for a different velocity, successive returns are most proximal. The term "most proximal" in the context of the present invention should be taken broadly and to mean "most likely to correspond to the expected signals reflected from a target traveling at that specific predetermined trajectory." These lines have different slopes in coordinates of range and time of transmission of the RF pulses, which coordinates are orthogonal to each other. An array of such lines is called a "butterfly pattern", and the lines are called "butterfly lines". The determination of which butterfly line the returns fall on or are most proximal to is called a "butterfly search". By "RF" is meant sonic (including ultrasonic) and electromagnetic radiation generally; ultrasonic radiation being presently preferred for medical diagnostic applications.

The invention, briefly described, estimates the velocity of a moving target, such as blood flow, from pulse-echoes received in response to a plurality of repetitive transmissions of propagating pulse signals. The echoes are translated into a plurality of time displaced outputs which are located in a space having at least two dimensions defined by orthogonal axes, one of which is the times of transmission of the RF pulse signals and the other of which is the time elapsed between transmission of these pulse signals and the reception of the echoes in response thereto. The system operates by detecting upon which of a plurality of butterfly lines in an array containing that plurality of lines, the outputs corresponding to the echoes are most proximal. The detection of the most proximal line may be carried out using a butterfly search. Each line has a different slope and defines a delay trajectory of the times of return of the echoes. The trajectory of the echoes (their time delay relationship) is effectively compared with the butterfly line trajectories. The estimate of the velocity of the target (e.g, the velocity of the ultrasonic energy scatterers which move in blood flow) is provided by the line having the trajectory most proximal to the trajectory of the output pulses.

The echoes may be translated into signals by processing the RF returns and utilizing the RF bursts. The envelopes of the RF bursts may also be used, or the bursts corresponding to the returns may be quadrature demodulated to provide quadrature related components as the outputs. These outputs are preferably digitized and sampled. Different signals constituted of samples which fall along the different butterfly lines are obtained. The butterfly line whose samples most closely correspond to the expected echoes returning from a target traveling at that predetermined velocity, is selected as the most proximal butterfly line. The predetermined velocity of this butterfly line is output as the velocity estimate. In the case of an envelope signal or an RF signal, a constant amplitude is expected from a butterfly line that exactly matches the trajectory of a target traveling at constant velocity. In the case of quadrature components, a constant amplitude, sinusoidally varying signal (of known frequency), is expected from a butterfly line that exactly corresponds to the echoes from a target traveling at constant velocity. In either case, all butterfly line signals can be evaluated to choose the single line with maximum correspondence between reflected echoes and the expected signal from a target traveling at velocities corresponding to the slope of each butterfly line.

The foregoing objects, features and advantages of the invention, as well as presently preferred embodiments thereof, will become more apparent from a reading of the following description in connection with the accompanying drawings in which:

FIG. 1 is a plot showing in slow time (pulse repetition period) and fast time (range) coordinates, an array of butterfly lines and ultrasonic RF return pulses which fall on one of the lines for a target velocity v.sub.0 ;

FIG. 2 is a plot similar to FIG. 1 showing the demodulated RF pulse envelopes of the returns;

FIG. 3 is a view similar to FIG. 2 but showing the butterfly lines in inverted "V" form and intersecting on the line corresponding to the first transmission;

FIG. 4 is a plot similar to FIG. 2, showing only one of the quadrature components and showing the samples obtained for echo returns, the maxima of which is shown by way of example as falling along one of the butterfly lines for velocity v.sub.0 ;

FIG. 5 is a block diagram of an ultrasonic imaging system having facilities for conducting a butterfly search to obtain the velocity of transcutaneous tissue, such as blood flow velocity;

FIG. 5a is a block diagram of a butterfly sampler shown in FIG. 5;

FIG. 6 is a block diagram of the maximum correspondence selector shown in FIG. 5 for the case where quadrature components are processed;

FIGS. 7a & b are respectively block diagrams of an embodiment of the maximum correspondence selector system of FIG. 5 using an envelope detection demodulator, sampler array and butterfly line variance detection array which is useful when RF pulse returns or the envelopes thereof are processed in a butterfly search for the most likely target velocity;

FIG. 8 is a block diagram of another embodiment of the system which produces quadrature components and uses complex processing to premultiply the complex envelope (quadrature I and Q components) with a complex operator, the sampler and a maximum correspondence selector of less complexity than the estimator of FIG. 6 and which may be used in the system shown in FIG. 5;

FIG. 9 is a block diagram of another embodiment of the system using complex processing of the RF returns with a Hilbert transform operator, the system being otherwise similar to that shown in FIG. 8;

FIG. 10 is a plot showing the variance detector outputs from of the RF envelope signals on the butterfly lines for an exemplary case of four transmissions (n=4) and with no additive noise which shows that for this case the variance minimizes on a butterfly line corresponding to the correct velocity of -75 cm/sec. (minus velocities indicating a target approaching the point of transmission of the ultrasonic pulses and positive velocities indicating a target moving away from the point of transmission); and

FIG. 11 is a plot showing the correspondence detector outputs for different butterfly lines for an exemplary case of n=4 returns and showing that the correspondence detector outputs with no additive noise, maximizes on a butterfly line for the correct velocity of -75 cm/sec. (minus velocities indicating a target approaching the point of transmission of the ultrasonic pulses and positive velocities indicating a target moving away from the point of transmission).

Referring to FIGS. 1, 2, 3 and 4, there are shown plots in two dimensional space along two different time axis which are orthogonal to each other. The abscissa is in much slower time than the ordinate. The abscissa represents successive Doppler pulse transmissions of ultrasonic, RF signals. The transmissions are separated in time by the reciprocal of the pulse repetition frequency which is the pulse repetition period or interval. The repetitive transmissions are n=0, 1, 2 . . . N. There are five transmissions in the plots of FIGS. 1-4. A feature of the system provided by the invention is that the velocity of the target is determinable with two transmissions. A practical and accurate system needs at least three transmissions. A set of eight or less transmissions (repetitive RF pulses transmitted towards the target) is presently preferred.

The ordinate axis is the fast time or twice the range from the location where the transmission is originated (the round trip travel time of the RF wave). The entire set of returns or echoes preferably occurs in less than the period between repetitions (PRT). Therefore, the ordinate is in fast time relative to the abscissa which is in slow time.

Drawn on the plots is an array or pattern of intersecting lines which are the butterfly pattern of butterfly lines. The intersection is at a particular range where the target is located; for example, the distance between the transducer which transmits and receives the ultrasonic, RF pulses and the blood vessel where the blood flow velocity is to be measured. The pattern may be effectively a "V" pattern, such as the inverted "V" shown in FIG. 3. There, the range is closer than in FIGS. 1, 2 and 4 to the target of interest. The location of the intersection and the butterfly pattern may be varied automatically, utilizing the range determining components of an ultrasonic imaging system embodying the velocity detection system described herein, as will be explained in connection with FIG. 5 hereinafter.

Each line of the pattern has a different slope and corresponds to a discrete linear velocity. The line to which the echoes are most closely proximal, and in the exemplary case as shown in FIGS. 1-4 where, for convenience of illustration, the peaks of the RF return or the demodulated versions thereof occur, represents the velocity estimate. This line is indicated as v.sub.0 for convenience. If all of the pulses were on the line indicated as v.sub.-2, the target reflector would be stationary since the travel time for each transmitted pulse and its echo then the same. In general, positive velocities away from the source (the transducer) are represented by butterfly lines of positive slope, while velocities of targets moving towards the source are represented by lines of negative slope as defined by the (n,t) coordinate system.

The correct velocity estimate is obtained by a butterfly search. This search is conducted to determine whether the time relationship of the echoes or returns is along the trajectory in the two dimensional, slow time--fast time space, defined by a butterfly line. As the description of systems implementing the butterfly search will make more apparent, the trajectories are sampling trajectories for echo samples in such time relationship that they fall on or are proximal to one of the lines in the butterfly pattern. Each sampling trajectory results in a one dimensional signal for all of the samples,

where i denotes each butterfly line in the array, that is the lines for v.sub.1, v.sub.0, and v.sub.-1 to v.sub.-5 in the plots illustrated in FIGS. 1 to 4.

FIG. 1 illustrates the echoes as bursts or pulses of ultrasonic RF signals, each being a US (RF) ultrasonic RF pulse. FIG. 1 shows the correct estimate is v.sub.0 since the expected constant amplitude samples of each of the five RF pulses 0 to 4 coincide with and fall on the butterfly line corresponding to v.sub.0.

FIG. 2 shows the use of the RF pulse envelope after demodulation. Just like the constant amplitude of the RF pulse shown in FIG. 1, the constant amplitude of the envelope falls on the v.sub.0 butterfly line. The constant amplitudes of both the RF pulse and the envelope also have minimum variance and, indeed, are constant in FIGS. 1 and 2 along the butterfly line v.sub.0. There is a constant delay .tau..sub.0. This delay is a constant travel time for the transmission and its echo, and may be subtracted from each round trip travel time to obtain the time delay relationship of the echoes. This time delay is, in the plots of FIGS. 1 to 4, on the trajectory of the v.sub.0 butterfly line. The butterfly search finds which of the butterfly lines corresponds to the trajectory of the echoes as defined by their time delay relationship. Thus, by matching the sampling trajectory of the butterfly lines with the actual delay trajectory of the echoes, the line that has the most constant amplitude or minimum variance will give the correct estimate of the velocity of the target, in effect the maximum correspondence to the butterfly lines is obtained. The processing is somewhat different for rf and envelope from that for quadrature components as will become more apparent hereinafter.

FIG. 4 illustrates one quadrature component resulting from quadrature demodulation of the RF pulse echoes. The combination of both quadrature components, along the butterfly line corresponding to target velocity v.sub.0 is defined as

both components of

vary in amplitude, but does so as a constant amplitude sinusoidal function. If the velocity of the target changes, so would the correct butterfly line and the frequency of the sinusoidal function. In the butterfly search using quadrature components, correspondence functions or ratios are obtained from each set of samples. Since they maximize on the correct butterfly line, they indicate the correct butterfly line, such that the velocity corresponding thereto, may be taken as the velocity estimate of the target.

In FIGS. 2 and 3 where envelopes from the signal reflector are shown, this search is along butterfly lines of constant velocity in two dimensional (fast time, slow time) space. As the reflecting target moves, the echoes are translated in time. Samples along a line of the correct velocity (the v.sub.0 butterfly line in FIGS. 2 and 3) have minimum variance. The butterfly search is carried out by sampling at different delays for each of the n echoes. Whichever butterfly line produces the minimum variance (the v.sub.0 line in the plots), is the correct estimate of the reflector target velocity. As pointed out above, in FIGS. 2 and 3, where the envelope has Gaussian variation, on the v.sub.0 line, all samples have the same value and hence, the variance is zero. In addition to minimum variance, the correct line may be found using minimum sum of absolute difference computations. A finite number of butterfly lines are used. The number of lines and their spacing may be increased or decreased depending upon the desired resolution in the velocity estimates. A large number of lines may be needed in some cases. For example, lines corresponding to velocities from--120 cm/sec to +120 cm/sec may be suitable for most diagnostic, blood flow measurement applications. See FIG. 10 for the output of the processing of the envelope components of FIGS. 2 & 3 and FIG. 11 for the output of the processing of the quadrature components as shown in FIG. 4.

Referring to FIG. 5, a transducer 10 emits an ultrasonic RF pulse, for example, from 500 KHz to 50 MHz, preferably 3 MHz for blood flow measurements. This pulse interacts with a reflector or scatterer from the moving target where a part of the wave is reflected. The reflection due to each successive transmission (a transducer firing or ping) take more and more time to reach the transducer. If the target is stationary, the time difference between firing and receiving the reflection (the echo) will remain constant. If the target is moving toward the transducer, this time difference decreases progressively on each firing. The ultrasonic RF pulses are generated in a transmitter (TRANSM.) 12. The output from the transducer is an electrical RF pulse signal which is applied to a receiver 14. Gating circuits (not shown) may be used to prevent the high power transmitted pulse from reaching the receiver. Timing pulses (tmg) and the pulse received by the receiver designate the travel time and are utilized in a range detector 16. The output of the range detector provides the coarse range upon deducting the .tau..sub.0 value and provides a range shift which is used to move the butterfly trajectories to the correct range for the target of interest. The range output may also be applied to a brightness or B display for displaying an image when the ultrasonic pulses are scanned over a target area, for example, the heart or womb, to provide an image of an organ or embryo of interest. To limit the search for velocities in ranges of interest, it may be desirable to change scale and filter settings which change the number of butterfly lines, the spacing between the butterfly lines and the range of search, and thereby shift the velocity ranges which are covered by the lines.

The echo signals from the receiver may be applied to a pre-processor which includes moving target indicator (MTI) rejection circuits thereby reducing clutter and noise in the signals, s[n,t] which are presented to the detection and velocity sampling and estimating elements of the system. A preprocessor 20 with conventional MTI or clutter elimination circuitry may be used for this purpose. The echo signals from a point spread target, have the form and can be expressed in terms of the first echo, in a set of echoes by resetting the time axis origin each time the transducer 10 fires. This echo signal may be expressed, in simplest form, neglecting attenuation, beamwidth, and other system effects, as ##EQU1## In discrete time representation for sample i, this expression becomes ##EQU2## In equations 2.01a and b: A is the signal amplitude (a constant); n is the index for the repeated echoes (the slow time index, n=0, 1 . . . 4 in the plot shown in FIGS. 1, 2 and 4); T is the pulse repetition period, .phi..sub.0 is a phase term depending upon the distance of the target from the transducer 10 at the time of first firing (n=o); .omega..sub.o is the angular center frequency of the RF ultrasonic transmission; r(t) is the envelope of the transmitted pulse, v is the velocity of the target which is v.sub.0 and can vary over the range from -v.sub.p to +v.sub.p (.+-.v.sub.p being the velocity range of interest); f.sub.f is the fast time sampling frequency where the i samples are taken, c is the velocity of the ultrasonic waves in the target medium, and the factor of two in equations 2.01a and b comes from the round trip travel time of the wave.

Three different types of signals may be processed and any may be selected. For the RF signal processing, where the RF pulse echoes shown in FIG. 1 is utilized, a direct connection is made. Otherwise, the signal may be demodulated by quadrature demodulator 22 or in a envelope detector 24. The resulting signals after envelope detection are denoted as e(n,t) and the complex envelope, consisting of the quadrature components, is r(n,t). The form of the RF ultrasonic pulse is the same as s[n,t] as given in equations 2.01a and b. For the envelope e(n,t) for the n lines ##EQU3## In discrete time representation ##EQU4## The signal of interest is sampled by the butterfly samplers for each of the butterfly lines. In FIG. 5 separate samplers 26, 28 and 30-32 are used for each butterfly line. Each sampler may be of the type shown in FIG. 5a.

For the ultrasonic RF pulse or the RF envelope, this search is to find the butterfly line over which the butterfly sampled signal is the most constant over increasing n. Then for the RF envelope, the quantity in the brackets in equation 3.01a should be constant. ##EQU5## In other words, the time difference between successive lines should be ##EQU6## In discrete time form ##EQU7## Since i can have only integer values, sample points in a butterfly line may lie between two successive values of i. In such cases, linear interpolation may be used to obtain the correct velocity estimate.

As shown in FIG. 5a, a multiplier 34 provides as an output, one dimensional signal for all the samples obtained around the range shift .tau..sub.o which are generated in a sampling pulse generator 36. The sampling pulse generator generates a delta function ##EQU8## In this delta function for the butterfly sampler 26, v is v.sub.-1. In the butterfly sampler 28, for the v.sub.0 line, v is v.sub.0. In the butterfly sampler 30 for the v.sub.-1 line, v is v.sub.-1. Thus, there are different sampling times which dictates separate searches along each of the butterfly lines in the butterfly array. The output of each sampler, when envelopes are processed, is a one dimensional signal e.sub.Bv [n] where v is the velocity for each of the lines. These one dimensional signals are applied to a maximum correspondence selector 38 to select the most likely velocity. Different estimators 38 are used depending upon the output signal which is sampled, and whether it is the RF echo, the RF envelope, or the quadrature components.

In the case of envelope signals, the maximum correspondence selector 38 selects the butterfly line with minimum variance and its slope can be used for the target velocity estimation. ##EQU9## where var(e.sub.Bv [n]) is the variance of the butterfly sampled envelope e.sub.Bv [n] on the v butterfly line.

When the output signal is the result of quadrature demodulation, the search for the correct velocity combines aspects of both time domain and frequency domain processing. The quadrature signals I and Q along the correct butterfly line will have a single frequency waveform whose frequency depends on the velocity, or slope of the butterfly line. Even in the presence of noise which would introduce other frequencies, the correspondence with the expected frequency will maximize on the correct butterfly line. The search is conducted where the I and Q components in each line is checked for the corresponding unique frequency as illustrated in FIG. 4.

As previously stated in equation (2.01a), the nth scan line will have the form ##EQU10## In discrete time representation per equation [2.01b] ##EQU11##

The envelope r(t) is assumed to be Gaussian, hence the quadrature components along any given line will be Gaussian. The complex envelope is ##EQU12## In discrete time representation, ##EQU13##

In case of the discrete time sampled data, for points lying between data samples, linear interpolation could be used.

Now, if we sample the complex envelope with the restriction, as imposed for the butterfly method on RF or envelope, that is, the time difference for the sampling trajectory between successive lines, ##EQU14## then, the complex envelope along this sampling trajectory is given by, ##EQU15## where .tau..sub.o is the reference time delay of the echo (FIG. 7), and r.sub.Bv [n] denotes sampling along the butterfly line for velocity v, ##EQU16##

As depicted in FIG. 4 (which illustrates the technique as applied for a single point spread target without noise), both of I and Q samples along the correct butterfly line would sample a single frequency constant amplitude sinusoid with respect to the index n. And, for each velocity, on the correct butterfly line, that distinct single frequency is unique, and under ideal conditions all the signal energy on the correct butterfly line would be concentrated only at that frequency. With noise present, that particular frequency would still be expected to have maximum energy. If we investigate each butterfly line for the power contained in the frequency corresponding to that line, normalized by the total power, the maxima should occur on the correct line (v=v.sub.O).

To illustrate this, consider the case of FIG. 4 and examine samples along other butterfly lines.

On an incorrect line, which would have sampling instants, ##EQU17##

As r() varies with n, these signals are now amplitude modulated and should be rejected.

We determine the frequency-specific energy along each butterfly line by applying orthogonality.

Since equation (3.08a) shows that the signal r.sub.Bvo [n] is a constant amplitude, single frequency complex exponential, and since the frequency term ##EQU18## is known a-priori and is unique for this trajectory, then by orthogonality, ##EQU19## Thus, let ##EQU20## Then ##EQU21## and, ##EQU22## where an incorrect choice of butterfly line (and velocity) will produce a lesser outcome in comparison to the total energy of

For the butterfly search over quadrature signals, a likelihood function L(v) is evaluated in maximum correspondence selector, shown in FIG. 6 for v=v.sub.0. To determine the likelihood, the samples along each butterfly line are evaluated by maximization thereof which occurs when the complex conjugate for the expected signal, if it would have occurred along that butterfly line, is multiplied with the actual signal. Accordingly, the complex conjugate is generated in a generator 40 and multiplied in a multiplier 42 with the signal. In FIG. 6, the sample signal along the v.sub.0 line is considered. Each echo of the n echoes provides a product which is summed in a summer 44. The sum is squared in a squaring circuit 46. The output is normalized by dividing by the reciprocal of the signal obtained from the v.sub.0 butterfly line sampler. That signal is squared in a squaring circuit 48. The sum of the squares is taken in a circuit 50 and inverted in a reciprocal generator 52. The product obtained in a multiplier 54 is the likelihood ratio or function, L(v.sub.0) of the signal on the v.sub.0 line.

All of the functions of the elements 40 to 54 may be performed by operations in a digital computer on the sample values or by dedicated hardware. The computer may be programmed in accordance with conventional programming techniques to provide the L(v) output. All the L(v) outputs are applied to a maximum detector 58 and the velocity v corresponding to the maximum L(v) will be read out as the correct velocity. This velocity may be applied to the display and used in imaging operations or in the spectral Doppler mode.

Mathematically, it can be shown that the maximum obtained from the likelihood functions correctly identifies the butterfly line corresponding to the correct velocity.

For each butterfly line we evaluate a likelihood term, ##EQU23## L(v) maximizes only at v=v.sub.0. So, we can determine the velocity by finding the maximum of L(v), ##EQU24## and, the mean velocity can be evaluated by finding the following weighted mean, ##EQU25##

The correct velocity may be obtained by premultiplication of the quadrature components with the complex exponent to simplify the likelihood ratio detectors for each butterfly line sample set. The system as shown in FIG. 8 uses such premultiplication. This simplification is obtained by evaluating the expression for L(v) given in equation 3.15 and repeated below.

For the estimation, we evaluate ##EQU26##

We can determine the velocity by finding the maximum of L(v), ##EQU27## and, the mean velocity can be evaluated by finding the following weighted mean, ##EQU28## Now, ##EQU29## Let us define,

Then, ##EQU30## which is just a summation of a sampled signal.

Thus, the likelihood ratio L(v) can be evaluated as follows, ##EQU31##

Simplifications in the likelihood ratio detectors which provide the L(v) outputs, may thus be obtained as shown in FIG. 8 by complex processing where quadrature demodulation is followed by fast-time multiplication by a complex function. Similar symbology is used as is used in FIG. 6. ##EQU32## From premultiplication, ##EQU33##

By comparing (4.13) and (4.14) we see that (4.14) can be obtained from the RF signal by just zeroing the positive frequency part in (4.13). This can further reduce the computation complexity, because this does away with complex envelope computation.

Thus, a similar likelihood ratio detector may be achieved even without the premultiplication step by using a Hilbert transform operator 60 as shown in FIG. 9. This step is identical to premultiplication by the complex exponential as shown in FIG. 8 and discussed above. That is why the steps after the complex processing are identical in FIG. 8 and FIG. 9.

FIGS. 7a and 7b illustrate the envelope detector 24 and butterfly samplers 26, 28 and 30 for the sake of simplification. The maximum correspondence selector in this case employs variance detectors 80, 82 & 84, which as shown for the detector 80 consists of a summer which sums all of the samples and a squaring circuit. The sum of the squares is evaluated by squaring the sampled signal e.sub.Bv [n] and then taking the sum of the squares. The square of the sums is subtracted from the sum of the squares in an adder circuit 62 to provide the variance signal. The minimum variance signal is obtained in a minimum variance detector 64 to identify the correct butterfly line. In other words, the minimum variance detector 64 is effectively a comparator or which selects the maximum output which defines the output velocity estimate v. A similar arrangement of variance detectors and comparators may be used to obtain the correct velocity butterfly line utilizing the RF ultrasonic pulses directly.

From the foregoing description, it will apparent that there has been provided an improved system (method and apparatus for velocity estimation) which is applicable to blood flow velocity measurement, as well as generally for measuring the velocity of moving reflectors. Color Doppler and color velocity imaging systems can utilize the improved system, whereby a butterfly search is conducted over each range of interest, and the velocity estimate for each depth is displayed as a color. Other variations, such as averaging of signals and averaging of butterfly estimates from nearby ranges in order to improve the noise immunity, are within the scope of this invention. Other variations, such as setting thresholds on the L(v) signals, below which the data are deemed unreliable may be applied. These and other variations and modifications in the herein described system, within the scope of the invention, will undoubtedly suggest themselves to those skilled in the art. Accordingly, the foregoing description should be taken as illustrative and not in a limiting sense.