Doppler Imaging Device

This application in part discloses and claims subject matter disclosed in our earlier filed pending application Ser. No. 07/891851, filed Jun. 1, 1992. This invention relates to acoustic imaging device, and more particularly, to Doppler imaging devices.

The prior art of Doppler detection and imaging techniques has advanced mainly in three areas: continuous wave Doppler detection, pulsed-Doppler detection, and color-flow Doppler imaging.

Continuous-wave Doppler systems continuously transmit ultrasound into human tissue. Echoes arise in the tissue and return to the receiving transducer. Dominating the return echoes are very strong signals from stationary sources. Superimposed on the strong signals are smaller Doppler shifted waves arising from moving tissue components.

CW Doppler systems normally use a single split element transducer, which continuously transmits and receives ultrasound. One transducer element transmits the carrier frequency into the tissue and the other receiving element amplifies the return signals, which are compared to the carrier signal to produce the Doppler signal. The Doppler signal then goes to a spectral analyzer which outputs the Doppler frequencies.

CW Doppler is used in medicine for monitoring blood flow velocities. The CW Doppler system obtains information along a line extending from the transducer. Thus, it does not produce an image of moving blood because of the lack of range information.

Pulsed Doppler imaging is a well known technique commonly used in ultrasonic medical applications. A transducer transmits a tone burst of acoustic energy into a medium and from the backscattered pulse interrogates the Doppler shift by mixing with the transmitted pulse. Range-gating, beam focusing and beam steering allow collection of Doppler information at a specific volume element inside the tissue. The resolution is the size of the volume element which is laterally determined by the focused beam waist and axially determined by range gating. Pulsed-Doppler tracks the phase of the received signal at the pulse repetition frequency (5-20 kHz), thus yielding a time series of the Doppler shift at a specific point.

Color Flow Doppler is pulsed Doppler displayed on multiple points to form an image. A linear phased array of transducers forms dynamically focused beams that scan the array. The image is separated into sampling sites approximately one acoustic wavelength long. The phase information is used to indicate motion, and based on phase changes the system separates an echo signal into forward and reverse channels, which depict the presence and direction of motion as colors.

A time-domain correlation method is also used to form an image of the velocity of moving blood. This method correlates in time successive pulse-echo scan lines to obtain an image of the velocity of moving blood. The direction and velocity of the moving blood is again depicted as colors.

Both color-flow Doppler imaging and time domain correlation imaging form the image first with a beamformer, and then look for changes in that image by using a Fourier transform or correlation technique.

The present invention provides a device and method for producing an image of a fluid moving through a medium. An acoustic wave insonifies the medium to produce Doppler shifted scattering from randomly located scatterers moving in the medium. At least eight channels of acoustic detectors detect signals reflected from the medium and the amplitude of the scattered signal is determined at each channel for at least one Doppler shifted frequency. Time-averaged receiver channel-to-channel amplitude correlations are produced using the Doppler shifted amplitude data from which the image of the moving fluid is calculated.

The invention disclosed in this patent application differs from the prior techniques in the following ways:

1. We form amplitude correlations between individual receiver channels.

2. We form these correlations without first using a beamformer to form an image.

3. We form an image from the receiver channel-to-channel amplitude correlations with algorithms not used in the prior techniques.

FIG. 1 is an illustration of the system of acoustic transducers utilized in this invention.

FIG. 2 is the spectrum of received Doppler information from simultaneous insonification at three discrete frequencies.

FIG. 3 is a schematic diagram of the invention used to test the concept in a water bath.

FIG. 4 is a schematic diagram of the invention used to image blood flow in a leg.

FIG. 5 is a detailed illustration of several of the acoustic transducers.

FIG. 6 is block diagram illustrating the electrical circuitry used in the processing of electrical signals from the transducers and the transmission of electrical signals to the transducers.

FIG. 7 is a block diagram showing the principal steps utilized in creating a Doppler image using a preferred embodiment of the present invention.

FIG. 8 is a simulated computer reconstructed image of flow in a cylindrical vessel using the mathematical processes according to the present invention. The simulation utilizes 128 transducers arranged regularly on a 120.degree. circular arc.

FIG. 9 is a simulated computer reconstructed image of flow in a cylindrical vessel using the mathematical processes according to the present invention. The simulation utilizes 64 transducers spaced regularly on a 120.degree. circular arc.

FIG. 10 is a computer reconstructed image using actual flow data. Image consists of flow through two tubes with 5 mm inner diameters and with equal but opposite flow velocities.

FIG. 11 is a computer reconstructed image using actual flow data. Image consists of flow through a single tube with 15 mm inner diameter.

One preferred embodiment of the present invention provides a system for insonifying a medium with a single frequency acoustic wave from a single direction, a system for measuring the amplitude and phase of acoustic waves scattered into a plurality of directions from the medium, and an efficient mathematical algorithm for determining a two-dimensional flow velocity map from this set of acoustic wave information. A schematic diagram of the cross section of a portion of this embodiment is shown in FIG. 1.

In this embodiment, there are N=256 acoustic transducers 10.1, . . . , 10.N arranged in the r-.theta. (or x-y) plane on the locus of a circular 360 degree arc of radius r.sub.o. The transducers 10 act as receivers of acoustic waves. The z-direction is perpendicular to this plane. The transducers circumscribe the medium 12 which is semi-transparent to the acoustic waves. In this preferred embodiment, a fluid resides in the remaining volume 11 to effectively couple acoustic waves between the transducers 10 and the medium 12. Water is a preferred coupling fluid for imaging human tissue. When imaging other objects a fluid should be chosen which preferably has a density .rho..sub.o and sound speed c.sub.o which are approximately equal to the average density and sound speed of the medium 12.

A continuous acoustic wave p.sub.o (t):cos[k.multidot.r-.omega..sub.o t+.phi.], where .phi. is a constant phase delay ##EQU1## and .omega..sub.o =kc.sub.o =2.pi. f.sub.o, where f.sub.o is typically 1 to several MHz, is transmitted from acoustic transducer 14 into human tissue as shown in FIG. 3. The tissue contains a plurality of vessels 13 which carry flowing fluid such as blood, for example. The moving blood cells scatter Doppler-shifted sound waves into all directions. The Doppler shift caused by a blood cell with velocity v is .omega.=k.multidot.v, corresponding to a frequency shift ##EQU2## For typical blood flow, v.ltoreq.1 m/sec giving f.ltoreq.1000 Hz for an acoustic wavelength .lambda.=1 mm (.function..sub.o =2 MHz) thus, the Doppler frequency f<<f.sub.o. The acoustic wave incident on receiver 10.j can be written p.sub.j (t)=m.sub.j (t)exp[i.omega..sub.o t], a form which separates the signal into its carrier component and the complex amplitude m.sub.j (t) containing the Doppler shift. The complex amplitude m.sub.j (t) contains a constant pan due to a combination of acoustic waves scattered from the static tissue and the time-dependent part due to Doppler-shifted acoustic waves scattered from the flowing blood.

The measured complex amplitude m.sub.j (t) at receiver 10.j can be written as:

where D.sub.j (r) is the antenna pattern for receiver 10.j. The source function .function.(r,t) is the complex amplitude for the transmitted acoustic wave to scatter from a volume element located at position r. The volume element can contain static tissue (static .function.(r)) or flowing fluid (time dependent f(r,t)). We Fourier transform equation (1) in time to yield

where F(r,.omega.) is the Fourier transform of the source function .function.(r,t). Next, we correlate all pairs of receivers by forming the correlation matrix

where <.....> denotes an ensemble average over different realizations of X.sub.j (.omega.). In the case of flowing fluid, a realization of X.sub.j (.omega.) consists of the acoustic Doppler spectrum from a specific collection of discrete scatterers. The algorithm relies on the fact that the scattering particles in the blood are randomly distributed in the blood vessel. Furthermore, since the particles are moving, the distribution of particles in the scattering volume varies randomly as a function of time, so a time average of the moving particles is equivalent to the ensemble average in equation 3. In the mathematical terms of random processes, the algorithm relies on the fact that the spatial part of the complex source function .function.(r,t) consists of a sum of localized point-like scattering points (with size <<1), each multiplied by a random complex amplitude exp[i.xi.], where .xi. is a Gaussian random phase. Thus, we can replace the ensemble average on the right-hand side of equation (3) with the form appropriate to random incoherent radiators:

In equation (4), S(r,.omega.) is the power radiated at r at Doppler frequency .omega.. Inserting equation (4) into equation (3), we arrive at

We will denote the ensemble average correlation matrix

Obtaining an image of the flow medium using our embodiment consists of inverting equation (5) for the source power distribution, S(r,.omega.).

The antenna pattern D.sub.j (r) appearing in equations (1), (2), (3), and (5) represents the response of the receiver 10.j to incident acoustic radiation. The radiation pattern of any acoustic transducer can be exactly calculated as a weighted sum of monopole, dipole, quadrupole, . . . patterns, as discussed, for example, by J. D. Jackson in Chapter 16 of "Classical Electrodynamics" (John Wiley, New York, 1962). Our Doppler reconstruction algorithm can be executed with any choice of D.sub.j (r), including the experimentally measured antenna pattern. For our preferred embodiment, the transducers closely resemble dipoles in the x-y plane (see FIG. 1), so D.sub.j (r) is the solution of the Helmholtz equation with dipole source term:

The antenna pattern in the z-direction is collimated in our preferred embodiment.

In equation (7), ##EQU3## corresponding to a coherent sum of two monopole (point) receivers which straddle the boundary along the radial direction. Using for the polar geometry of FIG. 1, the solution to equation (7) is ##EQU4##

In equation (8), we have omitted a normalization constant which does not effect the reconstruction algorithm, since only the distribution of the data, and not the absolute magnitude, is important. The function J.sub.m (z) is the Bessel function of integer order m and H'.sub.m (z) is the first derivative of the Hankel function of order m. These functions are tabulated in Chapter 9 of "The Pocketbook of Mathematical Functions" by M. Abramowitz and I. Stegun (Verlag Ham Deutsch, 1984). In this reference one can also find Graf's identify, which allows equation (8) to be simplified to ##EQU5##

Equation (5) theoretically justifies our principal claim, namely that detection using the embodiment pictured in the location, r, FIG. 1 provides incoherent Doppler-scattered data that suffices to characterize the source power distribution of flowing blood in tissue. We will now proceed to discuss several ways to solve equation (5) for the source function S(r,.omega.). S(r,.omega.) is power spectrum of acoustic signals scattered at a Doppler shifted frequency .omega. and position r and thus presents an "image" of the flowing fluid.

In theory, we could reconstruct, from equation (5), the source power distribution S(r,.omega.) as a continuous function of the Doppler shift frequency .omega. at each spatial point r. In practice, we will reconstruct S(r,.DELTA..omega..sub..beta.) for an average of frequency bands .DELTA..omega..sub..beta. and finite spatial pixels r.sub..varies..

To discuss the inversion techniques, first we adopt a notation which allows us to write the equations of Doppler imaging in matrix-vector form. Let Greek indices (.alpha., .beta., . . . ) refer to coordinates r in the x-y or r-.theta. plane, so a function G.sub.j (r) labeled by the jth receiver coordinate and evaluated at the pixel location r.sub..alpha. is denoted by G.sub.j.alpha.. Integrals over the coordinate r are then written as a sum over .alpha.; for instance ##EQU6## Also, we denote the average of Doppler frequency bands as .DELTA..omega..sub.62 . This is depicted graphically in FIG. 5. With this notation, we can rewrite equations (2), (4), (5) and (6): ##EQU7## In equation (10), we have added noise to each measurement X.sub.j, so the covariance of the noise, N.sub.ij.ident.n.sub.i n.sub.j, now appears in equation (12). In equations (10)-(12), it is implied that the variable .omega. is fixed at an average of a specific frequency band .DELTA..omega..sub.62 . We will reconstruct spatial maps of S(r.sub..varies., .DELTA..omega..sub..beta.) for each separate frequency band .DELTA..omega..sub..beta. and then combine this information to form an image.

The above equations define a class of mathematical problem called beamforming, which refers to the determination of a distribution of radiators given an ensemble of noisy field measurements along an antenna array. Equation (12) can be written in matrix form

A number of different techniques have been developed to solve this equation for passive beam forming applications, any of which could be applied to obtaining an image S(r,.omega.)=S.sub..alpha. (.omega.) of flowing blood from equation (13). We will briefly list some common methods, referring to the article "Multiple Emitter Location and Signal Parameter Estimation," by R. O. Schmidt (IEEE Trans. Antennas and Propagation, AP 14. p. 276, 1986) for comments on their relative merits. Then we will elaborate on one method, maximum likelihood, which we have shown to work well on real data. Finally, we will also show that a reconstruction method that we developed for coherent ultrasound tomography (see pending patent application Ser. No. 07/891,851, filed Jun. 1, 1992) can be used to obtain a Doppler image.

Ordinary beamforming performs a direct inversion of equation (13):

where the subscript .alpha..alpha. denotes the diagonal term of .alpha.column of the square matrix.

Because the antenna pattern is broad, S.sub..alpha..sup.BF will show broad peaks in the direction of blood scatterers, but this method does not resolve different scatterers as required for an image. The so-called music algorithm (referenced above by Schmidt) solves equation (13) for the case of sparse radiators by separating the vector space into a signal space (spanned by the vectors D.sub.j) and a noise subspace (spanned by the eigenvectors of N). The number of incident signals can be determined by computing the eigenvalues of R in the noise representation, and the image is then given by ##EQU8## where E is the matrix whose columns are the eigenvectors of the noise matrix, N. In non-medical applications, the music algorithm has been shown to locate distinct acoustic radiators in the presence of noise. The maximum entropy estimator obtains an image for each receiver, ##EQU9## where C is a column vector of R. It selects one of the array elements 10.j as a reference and finds weights to apply to the remaining signals to minimize the mean-square error between the reference signal and the linear estimate. The maximum likelihood estimator assumes that the measurements X in equation (13) correspond to a zero mean, multiparameter Gaussian process. Baye's theorem is used to derive an expression for the conditional probability, P(S;X), of a given source S given the measurements, X. The function P(S;X) is then maximized by minimizing the log-likelihood function, ln P(S;X), with respect to S, with the constraint Set.sub..alpha. >0, for all .alpha.. Conventional maximum likelihood uses ##EQU10## to obtain an image. We have found that equation (17) produces good images using our embodiment. Improved images can be generated by including higher-order correlations through extending the maximum likelihood formalism. Minimizing ln P(S;X) results in the following non-linear equation to be solved for S:

Equation (18) can be solved iteratively by assuming an initial first-order solution XX.sup.T =R.sub.1 =N+DSD.sup.T, so that (N+DSD.sup.T).sup.-1 .perspectiveto.R.sub.1.sup.-1. The first-order estimate for S.sub..alpha. is then

This solution has been shown to work well on our simulated Doppler data. An alternative iterative method follows if one uses (N+DSD.sup.T).sup.-1 =R.sub.1.sup.-1 only once in equation (18):

Equation (20) may be iterated by using the maximum likelihood estimate equation (17) (first term on right-hand side of (20)) as a first approximation S.sub.1.sup.-1 to , and S.sup.-1, and evaluating the second term on right-hand side of equation (20) using S.sub.1.sup.-1.

Another method of inverting equation (13) exploits the mathematical similarity to the reconstruction of coherent ultrasound tomography as described in pended patent application Ser. No. 07/891,851, filed Jun. 1, 1992. Ultrasound tomography solves for a source function S(r) given measurements

where j and k are transducers 10.j. The form in equation (21) results from the fact that in coherent tomography, receivers 10.j are used one at a time as transmitters, allowing measurements (21) to be mapped out for all pairs of transmitters/receivers. In equation (21), .alpha..sub.j (r) is the solution to the Hemholtz equation for a dipole transmitter at r.sub.j, including effects of multiple reflections. In the absence of multiple reflections, it is identical to the antenna function D.sub.j (r). Formally, equation (21) is identical to the Doppler equations (5)-(6), with M.sub.jk .fwdarw.R.sub.jk, .alpha..sub.j (r).fwdarw.D.sub.j (r), and .alpha..sub.k (r).fwdarw.D.sub.k *(r). Thus, the linear Fourier transform method described in the above patent application to solve equation (21) can directly applied to obtain an image from Doppler-scattered waves using our embodiment.

The frequency variable .omega. in the equations above corresponds to the Doppler frequency shift imparted to the signal scattered from that portion of medium 12 which is moving. The magnitude of the .omega. depends on the vector velocity of the medium, the vector wavenumber k, and the position at which the scattered wave is measured. For the case where the scattering medium is measured in a single plane perpendicular to the velocity of the medium .omega.=V.multidot.k. Thus, the solution, S(r.omega.), for equations (5) and (6) represents a statistical estimate of the expected or average number of particles at r with velocity corresponding to a Doppler frequency shift .omega.. In principle, a velocity map of all points in the medium may be determined by calculating S(r,.omega.) for all reasonable .omega.. In this case the following relation holds, ##EQU11## where V.sub.z is the average value of the velocity perpendicular to the plane containing the receiver transducers, C is a constant set by the configuration of the test device, and .theta. is the angle between the fluid velocity vector and the wavenumber vector k.

A means of improving the accuracy and reliability of the estimate V.sub.z (r) can be achieved by measuring R.sub.jk (.omega.) for more than one insonifying frequency .omega..sub.o. This could be accomplished sequentially; first insonifying the medium with a single frequency .omega..sub.o.sup.(1), measuring the scattered signal and storing the data, then repeating this step for as many frequencies, .omega..sub.o.sup.(i), as desired. Equations (5) and (6) could then be solved for each insonifying frequency and the final result would be the obtained by averaging the individual estimates of S for frequency shifts corresponding to the same scattering particle. If the total number of insonifying frequencies used is I then the velocity may be written as follows: ##EQU12## In practice, using multiple insonifying frequencies will help alleviate errors arising from the inexact knowledge of the antenna functions D.sub.j (r) .

If the velocities of the scattering medium are confined to a finite range such that .parallel.v.sub.m .parallel.<<c, where c is the soundspeed of the insonifying wavefront, and v.sub.m is the maximum velocity present in the flow, then the Doppler shifts are relatively small compared to .omega..sub.o. In this case S(r.omega.) is narrowband, ##EQU13## and it is possible to transmit multiple insonifying frequencies simultaneously. The medium is insonified with a complex waveform which has power at several .omega..sub.o.sup.(i) separated in frequency by more than the largest possible Doppler shift. The scattered signal may then be separated into distinct wavebands; one for each of the different insonifying frequencies. For example, if the largest velocity is v.sub.m =3 m/sec, then the maximum Doppler shift is ##EQU14## at insonifying frequency ##EQU15## Therefore, one could simultaneously transmit a comb of frequencies, .omega..sub.o.sup.(i) separated by ##EQU16## FIG. 2 illustrates this situation for an insonifying signal composed of three sine waves. The power spectrum of the scattered signal at a single transducer 10.j is shown. There is a strong spike at each insonifying frequency due to the sound scattered from the stationary part of the medium 12. The side lobes to either side of the main insonifying frequency peaks represents the Doppler shifted waves scattering from the moving portions of 12. Since the .omega..sub.o.sup.(i) are separated by more than the largest Doppler shift, the side lobes around any of the main frequency peaks goes to zero before the side lobes from any of the other insonifying frequencies become nonzero. The determination of the scattering potential and fluid velocity can then be accomplished using equation (23) just as if the medium had been sequentially insonified with simple single frequency wavefronts.

A sketch of an embodiment of the device used to test the Doppler imaging concept is shown in FIG. 3. A continuous-wave 500 kHz sound wave is transmitted into a medium containing plurality of vessels which carry flowing fluid such as blood, for example. Inhomogeneities in the flowing fluid (such as red blood cells) scatter the sound wave in all directions. A separate receiver array 10 interrogates the scattered sound waves emanating from the fluid. In the present system the receiver array consists of 256 transducer elements evenly arranged on a circle. Our test set up used a water bath to couple the sound waves to a Latex tube 15 containing the moving fluid. The receiver array was connected electronically to a computer 17 for high-speed data acquisition and data processing, and an image display unit 19.

The transducer portion of the present embodiment of the invention is displayed in FIG. 5. The device consists of N=256 acoustic transducers 10 evenly spaced on a circle of radius r.sub.o =51 mm which act as receivers, and an insonifying transducer 14. The insonifying transducer is placed above or below the receiver transducers 10, and positioned in a way to insonify the test medium with a beam of ultrasonic radiation. Each of the acoustic transducers 10 performs as a receiver of the acoustic waves scattered by inhomogeneities with in the test medium. The frequency of the insonifying waves is set at approximately 500 kHz which results in a wavelength .lambda..sub.o .apprxeq.3.0 mm in water and nearly the same wavelength in human tissue. The spacing between the receiver transducers 10 is 1.3 min. In the preferred embodiment, the acoustic antenna pattern for each transducer is fairly uniform with angle but is somewhat more sensitive in the radial (r) dimension which lends itself either to the monopole or to the dipole approximation for the steering functions D.sub.j used in the aforementioned derivation of the reconstruction algorithm. The vertical height of the transducers is 25.4 mm, which produces a fairly collimated receiver antenna pattern in the z-direction. The thickness of the image slices is roughly equal to the vertical (z) height of the transducers so the image S(x) will comprise an average over the vertical dimension. Alternate embodiments include varying the vertical height to obtain different image slice dimensions.

Each transducer element 10.j (j=1, . . . , 256) is constructed from PZT5A piezoelectric ceramic supplied by Vernitron Piezoelectric Division, Bedford. Ohio. The dimensions of each element are 0.90 mm (3-axis).times.2.80 mm (1-axis).times.25.4 mm (2-axis). Electrical energy is supplied along the 3-axis and the 1-axis is dimensioned so that the transducer can resonantly receive acoustic energy in a range centered around 50) kHz along the 1-axis. Each transducer backing 31 and the spacing between the transducers 32 is comprised of an air-filled silicone foam (Silpak SF-2000) which acoustically and electrically isolates the transducers from each other. FIG. 5 shows that one electrode of each transducer is connected to a common electrical ground 35 and the other electrode of each transducer 10.j is connected to an electrical strip line 42.j which is etched on a printed circuit board 40. Electrical signals are routed to and from each transducer along these strip lines 42.

A block diagram of the electrical processing circuitry is illustrated in FIG. 6. The electrical lines from the 256 transducers 10.j are connected through the electrical strip lines 42 to the inputs of sixteen 16-to-1 analog multiplexers 60.1, . . . , 60.16 where they are multiplexed. Each of the 16 output lines of these multiplexers are again multiplexed to a single line via another 16-to-1 analog multiplexer/amplifier 65. The output of this multiplexer 65 is sent to an analog-to-digital converter 68 which is synchronized with the signal from the signal generator 50. The digitized data from the analog-to-digital converter is stored in a computer 72 which for our prototype was a conventional 486 personal computer.

This invention can produce a two-dimensional map of the magnitude of the scattering potential of a medium moving within tissue which is semi-transparent to acoustic waves. The operation of this invention can be divided into two steps:

1. Collection of scattered wave information in the form of a complex time series.

2. Numerical calculation of the magnitude of the scattering power S(r,.omega.).

The two operational steps are described in detail below.

As described above the raw data collected by the piezo-electric receivers are instantaneous real scalar measurements of pressure at each transducer. In theory, a simple set of time series of the pressure at each receiver 10.j can be used to compute X.sub.j (.omega.), defined in equation (2) and the covariance matrices R.sub.ij (.omega.), equation (6), for any Doppler frequency shift. In practice however, the amount of data, electronics, and computations required to resolve fluid velocities typical in a blood vessel from a set of N distinct time series makes this approach unreasonable. So, a priori knowledge about the narrowband nature of the sonic signal received at the transducers is used to reduce the amount of data required. For blood flow, .vertline.v.vertline..ltoreq.3 m/s while c.apprxeq.1500 m/s, and since the Doppler frequency shift for a fluid traveling with velocity v is .DELTA..function.=.function..sub.o (v/c)cos.theta.. Therefore we know that the acoustic signal received at the transducers is narrowband with a bandwidth less than 2.function..sub.o v/c. And since p.sub.j (t)=m.sub.j (t)exp[i.omega..sub.o t ]

where the capital letter signify Fourier transforms. But, as described above, p.sub.j (t) is narrowband and therefore we need only sample m.sub.j (t) at a sample rate sufficient to accurately resolve a bandwidth of at most 2.DELTA.f.

The Doppler spectrum X.sub.j (.omega.) (j=1, . . . , N) for each transducer is collected by the embodiment in the following fashion. A single frequency sinusoidal electrical signal is sent from the signal generator 50 to the insonifying transducer 11 which transmits this signal as an acoustic wave into the medium 12. The transmitted acoustic wave scatters from the inhomogeneities of medium 12. The sound waves scattered from inhomogeneities which are moving relative to the insonifying transducer are therefore shifted in frequency. The scattered wave from the moving (blood) and stationary (other tissue) inhomogeneities impinge upon the N transducers 10.1, . . . , 10.256 which transform these acoustic signals into narrowband electrical signals. The electrical signals from the receiver transducers are sequentially routed to the analog-to-digital (A/D) converter 68 as shown on FIG. 6. The multiplexers are programmed so that the A/D converter will sample each transducer 10.j four times before the next signal, from the 10.(j+1) transducer, is routed to it. The sample rate of the A/D is set to ##EQU17## and the four samples are used to compute m.sub.j (t) for each channel. The result is stored in the computer's memory. This process is then repeated L times so that a time series of m.sub.j,l =m.sub.j (l.DELTA.T) is computed for each transducer, where (l=1,2, . . . , L) and .DELTA.T is the time elapsed between gathering the four samples required to compute m.sub.j,l and the four samples required to compute m.sub.j,l +1. For N=256 transducers

where T.sub.d is a user selectable delay time which can be varied depending on the magnitude of the anticipated Doppler shift, and t.sub.d is a fixed delay time which allows the electrical signal to settle after the multiplexer switches between channels.

A fast-fourier-transform (FFT) of length L is then computed on the data for each transducer yielding the quantities X'.sub.j (.omega.). The quantity X.sub.j (.omega.) is formed by multiplying the result of the Fourier transform by a complex value to compensate for the temporal offset between the time series for different receiver channels.

The result may be expressed as

This step is required because the individual transducers are sampled sequentially rather then simultaneously so that each element in the time series for the 10.j transducer is delayed by j(4.delta.t+t.sub.d). The covariance matrix R.sub.ij is formed for a particular frequency, shift .omega. by averaging the product X.sub.i (.omega.)X.sub.j *(107 ) for several sets of measurements--for the tests conducted with the prototype system 40-50 sets of measurements were averaged when generating the coefficients of R.sub.ij. Finally, the two-dimensional image of S(r,.omega.) is then generated using one of the reconstruction algorithms discussed previously.

Our computer software can be described by FIG. 7 which outlines the steps required to process the wave data for the image.

Step 1: Normalizes receiver data.

Step 2: Fast-Fourier-transforms m.sub.j,l

Step 3: Compensates for time delays between channels.

Step 4: Computes R.sub.ij from several sets of measurements.

Step 5: Reconstructs 2-dimensional image from R.sub.ij.

A computer simulation of the system was developed to test the reconstruction technique. This simulation uses an assumed antenna pattern D.sub.j (monopole, dipole or higher order function) to compute the signal at each receiver 10.j from a set or scattering points in the medium. Random noise commensurate to the noise expected in the real system is added to each measurement X.sub.j (.omega.). Each scattering point is assigned a fixed magnitude scattering power S(r,.omega.), and imparts a random phase shift to the scattered wave. Simulated measurements from several realizations of the random phase shift associated with S(r,.omega.) are used to form the covariance matrix R.sub.ij for a given Doppler shift. The reconstruction algorithms compute the scattering power S(r,.omega.) which is then compared to the known power distribution.

The results of the simulation are shown in FIGS. 8 and 9. FIG. 8 shows a reconstruction at a single velocity of simulated blood flow in a cylindrical vessel. A 120.degree. circular receiver arc is simulated with N=256. The flow field is assumed to be fully developed and laminar, with a velocity profile described by ##EQU18## where v.sub.o is the maximum velocity of the blood at the center of the vessel, a.sub.o is the radius of the vessel, and a is the radial dimension measured from the center of the vessel. The image is an annulus corresponding to fluid with uniform velocity at a radius a=.lambda..sub.o. FIG. 9 shows the same image reconstructed using only 64 receiver transducers arranged at uniform intervals in an arc of 120.degree..

FIGS. 10 and 11 show actual images obtained using the prototype device described above by reference to FIG. 3. The data was obtained with a full receiver ring of 256 receiver transducers 10. The transducer ring was mounted on a water reservoir such that the z axis (the normal to the image plane) was vertical. The reservoir was filled so that the transducer elements were submerged. The insonifying transducer 14 was placed above the receiver ring and was also submerged in the water bath.

FIG. 10 is a reconstruction of flow within the latex tube of inner diameter 0.5 cm shown in FIG. 3. The axes of the tubes 15 were perpendicular to the image plane and the flow within the tubes was equal in magnitude but opposite in direction: one flowing downward through the plane the other upwards through the imaging plane. Since the sound waves scattered by the fluid in the different tubes have different Doppler shifts, the two Doppler frequency bands were reconstructed separately to form the entire image. First the covariance matrix corresponding to a positive Doppler shift, R.sub.ij (.omega.)=X.sub.i (.omega.)X.sub.j *(.omega.), was used to reconstruct the image of the flow within the tube carrying the liquid upwards out of the reservoir. Next the data corresponding to a negative Doppler shift R.sub.ij (-.omega.) was used to reconstruct the flow within the other tube. To form the final image the two reconstructions were combined and the results are indicated on monitor 19 in FIG. 3. FIG. 10 is the actual computer generated image. Flow differences in color are shown on monitor 19.

FIG. 11 is a reconstruction of a single latex tube with diameter 1.5 cm. The flow within this tube is downwards into the reservoir.

Three-dimensional images of moving fluid can be obtained by reconstructing S(r,.omega.) as two-dimensional cross-sectional slices juxtapositioned in the z-direction and then stacking these slices to tom the three-dimensional image. The acoustic signal data used to reconstruct the two-dimensional slice images can be collected using a single one-dimensional array of transducers which is mechanically moved in the z-direction.

Alternatively, the acoustic signal data can be collected by using rows of one-dimensional transducer arrays arranged in the z-direction (i.e. a two-dimensional array). Again, acoustic data from each row of one-dimensional arrays is used to reconstruct S(r,.omega.) as separate two-dimensional cross-sectional images and then these images are stacked to form a three-dimensional image. Beam forming in the z-direction (i.e. phase weighted summing of individual transducer array elements spaced in the z-direction) can be used to reduce the slice thickness and increase resolution in the z-direction.

Finally, acoustic signal data can be collected with the entire two-dimensional array and this entire data set can be used to reconstruct S(r,.omega.) directly as a three-dimensional image.

Alternate embodiments of the invention also include the geometry of the transducer array. These alternate embodiments may be appropriate to facilitate coupling of the acoustic waves to various blood vessels in different parts of the human body including, for example, arteries and veins in the legs, arms, and neck.

The number and arrangement of the receiver transducers may be changed to fit the purpose of the device. The receiver array may be fabricated in either a circular ring as described above or a circular arc. Preferred design would provide arcs of at least 60.degree.. The number of receiver transducers may be varied from as few as 8 to more than 3072. The arrangement of transducers on the ring may be at regular intervals or random intervals. Designs which are not circular may be used, but the algorithms for reconstructing the image will be more complicated. FIGS. 4A and 4B depict an embodiment of the present invention using transducers arranged in a 120.degree. arc. This embodiment uses a fluid-filled bladder instead of a water bath to couple the acoustic signals into and out of the tissue. Persons skilled in the art of imaging devices can derive such modified algorithms using the principles described above.

An alternate embodiment concerns the vertical (z coordinate) dimension of the receiver transducers. The vertical dimension may be made smaller to image a narrower slice of the body or larger to average a larger slice.

An alternate embodiment concerns the nature of the insonifying transducer. The beam pattern emanating from this transducer may be broad to illuminate all the tissue inside the arc, or may be very narrow so that the operator may chose to look at only a select portion of the tissue. Another preferred embodiment would utilize a phased array of transducers so as to direct the insonifying beam to particular locations in the tissue.

An alternate embodiment concerns transmitting waves containing, multiple frequency components and using the added information to enhance the reconstructed images. Time-gating the reception of the multiple frequencyt pulses allows further discrimination of flowing fluid versus static tissue signals thus increasing the signal-to-noise ratio of the device. The reader should construe the above described embodiments of this invention as examples, and the scope of this invention shall be determined by the appended claims and their legal equivalents.