1022
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. IO. OCTOBER 1992
Calculation of Doppler Spectral Power Density Functions Geoffrey K. Aldis and Rosemary S . Thompson
Abstract-A volume integral method for the calculation of Doppler ultrasound spectral power density (spd) functions is described. Axisymmetric flow in a circular tube with a power law velocity profile is assumed. The spd function is regarded as a probability density function for scatterer velocity, and the assumptions under which this is justified are considered. It is shown that the spd function is independent of Doppler angle except in the presence of wall reflection effects. A coordinate system centered on the beam is used and this enables the integrals to be easily formulated for arbitrary beams. Irregularly shaped and nonuniform beams can be treated. For the common flow and beam patterns, which exhibit symmetry, the volume integrals can often be reduced to a single integral and evaluated directly. The method is applied and the spectra are calculated for various different cases. Results are obtained for uniform rectangular and circular insonating beams, and for nonuniform beams with Gaussian, jinc, and sinc profiles. The effects of narrow beams and wall reflection are shown. The method may be readily applied to other beam and flow patterns, and extension to more complicated situations is also discussed.
INTRODUCTION OPPLER spectra are readily obtained from different arteries in the body, but care is required in making inferences about physiological flow patterns from these spectra. The Doppler spectrum results from the flow velocity profile in the artery, but it may be modified by several other processes which are not flow-related. The spectral power density (the fraction of backscattered power contributed by scatterers at each velocity) may not be truly representative of the velocity density (the fraction of scatterers moving with each velocity). Following Evans et al. [l] the reasons for this can be summarized into three groups: i) scattering properties of blood or other target medium; ii) processing of the received signal by, for example, filtering or the FFT algorithm; iii) three-dimensional (3-D) geometry of the transduction process. In this paper the third factor, spectral distortion due to the 3-D geometry of the transduction process, is considered for continuous wave (CW) Doppler spectra. A num-
D
Manuscript received April 3, 1991; revised Februaty 10, 1992. This work was supported by the National Health and Medical Research Council of Australia. G. K. Aldis is with the Department of Mathematics, University College, University of New South Wales, Australian Defence Force Academy, Campbell ACT 2600, Australia. R. S . Thompson is with the Department of Obstetrics and Gynaecology, University of Sydney at Westmead Hospital, Westmead, NSW 2145, Australia. IEEE Log Number 9202575.
ber of different authors have previously discussed various aspects of this spectral distortion. Evans [2] and Evans et al. [ 11 considered the problem two-dimensionally using rectangular beams which may be narrower than the tube. Gill [3] examined errors in the estimated mean velocity, and presented some numerical calculations of spectra using a nonuniform beam profile. Cobbold et aZ. [4] developed a computer model of the 3-D problem incorporating effects due to a nonuniform beam profile, off-axis alignment of the beam, and differential attenuation between blood and tissue. Numerical results were obtained for various cases of interest and this approach was continued by Bascom and Cobbold [ 5 ] , who also derived an analytic expression for the spectral power density with a Gaussian beam profile. Cobbold’s model was also extended to study intrinsic spectral broadening by Bascom et al. [6]. Intrinsic spectral broadening is broadening of the received spectrum due to properties of the measurement system itself, and has been discussed extensively by other authors [l]. Newhouse et al. [7] and Censor et al. [8] presented theoretical results using the ultrasonic fields appropriate to certain transducers, whereas Bascom et al. [6] used a geometrical ray tracing approach. Spectral broadening can be described either geometrically or in terms of the finite transit time of scatterers through the beam [l]. Although it is a consequence of the 3-D geometry of the transduction process, spectral broadening is not considered here. This is equivalent to the assumption that all ultrasonic wavefronts are planar. Spectral broadening tends to smooth the shape of Doppler spectra [l] but since all frequencies are affected, this may not be evident except towards the extremities of the spectrum. The neglect of spectral broadening in this paper means that the calculated spectra will have an unphysical sharp cut off at the maximum frequency. In this paper we assume that the flow being insonated is steady axisymmetric flow of a simple power law form. We also assume that over the insonated region of the tube the beam has plane wavefronts and a fixed but not necessarily uniform cross-sectional intensity pattern. A circular transducer produces a beam for which the directivity function in the far field (Fraunhofer zone) is a jinc function, whereas for a finite length line source the directivity function is a sinc function [9]. The beam intensity pattern of ultrasonic transducers may be approximated using these
00 18-9294/92$03.OO 0 1992 IEEE
ALDIS AND THOMPSON: CALCULATION OF DOPPLER DENSITY FUNCTIONS
functions. A Gaussian beam intensity is also commonly used as it is simpler, and closely approximates both these functions in the axial region. Spectral power density (spd) functions can be calculated using volume fractions (Roevros [lo], Cobbold et al. [4]).A volume integral approach is developed in this paper using the 3-D geometry of the beam/tube configuration and an arbitrary beam profile. Integrals are formulated in beam coordinates which enables arbitrary beam profiles and shapes to be treated. It is shown that the spd function is independent of Doppler angle except in the presence of wall reflection effects. Even when the spd functions cannot be obtained analytically the volume integrals can often be simplified and computation reduced. A number of results are obtained using the method. Spectra are calculated for nonuniform Gaussian, jinc, and sinc beam profiles as well as for uniform rectangular and uniform circular beams. The effects of narrow beams and wall reflection are shown. METHODOF CALCULATION OF SPECTRAL POWER DENSITYFUNCTIONS
Assumptions The backscattered power from a small control volume of fluid is assumed proportional to the total number of scatterers in the volume and to the insonating beam intensity. We assume that there are a large number of independent point scatterers uniformly distributed throughout the insonated volume. Backscattered power is therefore proportional to volume if the insonating beam is uniform. When the beam has a variation in intensity over its crosssection the backscattered power is proportional to the volume weighted by the beam intensity. Over the insonated region of the tube the ultrasonic wavefronts are assumed perpendicular to the transducer face (no spectral broadening). There is no variation in the beam intensity pattern with axial distance from the transducer (thin-walled tube, no attenuation). We consider axisymmetric flow in a circular tube of radius a of the form
where q is the normalized velocity, r is radial distance from the tube axis, and n is a parameter. Poiseuille flow corresponds to n = 2 and as n -, 00 the profile approaches plug flow. All particles moving with velocity greater than q lie inside a flow cylinder of radius r = a ( l - q ) ' / " .The range of the normalized velocity q is 0 to 1 , but extension to any velocity range is easily accomplished by multiplying q by qmax.The Doppler spectra are presented in this study as functions of normalized velocity. The spectra could be equivalently regarded as functions of normalized Doppler frequency shift, since velocity and frequency shift are directly proportional by the Doppler equation.
1023
Spectral Power Distribution Function, F(q) Under the ideal conditions of uniform scatterer distribution and uniform insonation assumed here, the spd function is proportional to the velocity density and can be regarded as a probability density function for both spectral power and velocity. Under nonideal conditions the spd function may not be a probability density function for the velocity. Interest here is primarily on comparing the shape of observed spectra and not with the value of total backscattered power. The area under the spd function is therefore always left equal to 1 even under conditions of incomplete insonation. Associated with the spd function there i s a spectral probability distribution function, F(q) = j i f ( q ) dq. The value of the spectral power distribution function, F(q), is the fraction of the total insonated volume occupied by particles with velocities between 0 and q. F is a monotonic increasing function with F(0) = 0 and F( 1) = 1. It is often easier to determine F(q) and then to calculate the spd function from it by differentiation. The mean velocity, 4, can be obtained from Z j = Sh ( 1 - F(q)) dq = q f ( q ) dq. The median velocity qmis found from solving F(q,) = ( V I . The angle between the beam axis and the tube axis is the Doppler angle For plane wavefronts and axisym# metric flow a Doppler shift is only obtained for (a/2). However, it greatly simplifies the volume calculations if the transducer is inclined at a / 2 to the flow axis. We show below that for axisymmetric flow in a circular tube choosing Bo = (a/2) leads to the correct spd function. This holds even in the case of wall reflection, where the spd function depends on OD, provided the beam shape is appropriately modified. Calculation of F(q) Using Volume Integrals Integrals will be formulated using coordinate systems centered on the beam. The beam's cross-section has the general boundary A = A(%)and a uniform beam intensity over the beam cross-section is assumed for the moment. Fig. l(a) shows the beam and tube and their respective coordinate systems which share a common origin at the center of the tube. The tube Cartesian system (x', y', z ' ) is centered in the tube with the z '-axis along the tube axis and the y'-axis also in the plane of the page. The beam Cartesian system (x, y, z ) has its z-axis along the beam. The x- and x'-axes coincide. The y'- and z'-axes are obtained by rotating the y- and z-axes around the x-axis through OD. A cylindrical beam system ( r , 9, z) is also defined with x = r cos 9 and y = r sin 9. The surface of the tube is given by x f 2 y'* = a2 and since x ' = x and y ' = y cos - z sin 80 we find that the intersection of the beam and the tube occurs at
+
Fig. l(b) shows the projection of the tube onto the plane = 0 (a cross-section of the beam). The projection shows
z
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 39, NO. IO, OCTOBER 1992
1024
Hence
ir i
2
VT = sin 80 2'
rmax
r(a2 - r2 cos2
dr d€J ( 2 )
in which the dependence on Doppler angle is conveniently outside the integrals. Particles with velocities greater than q lie inside the flow cylinder of radius r = a ( l - q)'/".The volume occupied by the intersection of the beam with these particles can be obtained readily from the expression above for VTas
jr 1 r(a2(1
2 Vq = sin 60
'4
-
q)2/"- r2 cos2
dr de
(3) where
\
x=x'
X
(b) Fig. 1. (a) Cartesian coordinate systems for the beam ( x , y. z) and the tube (x', y', z ' ) . The tube system is obtained by rotating the beam system through Bo around the common x = 1'-axis. (b) A cross-section of the beam (in the plane z = 0) with the projection o f the tube shown shaded. Distance from the origin to the beam edge is A ( @ , distance to the projected tube edge is a / ( c o s 01.
the region of r and 0 which is used to calculate the total volume of intersection of the beam and tube. At each value of 0 integration extends in the r direction to either the edge of the tube, a l l c o s 81, or to the edge of the beam, A ( @ , whichever is encountered first. The total volume of the intersection of the beam and the tube, V,, is therefore
where
- (a2 -
r2 cos2
+ ( a 2 - r2 cos2
The volume fraction occupied by particles with velocities greater than q is V q / V p Therefore F(4), the volume fraction occupied by particles with velocities less than q , is F(q) = 1 - V q / V T . The dependence on Doppler angle cancels in F(q) for any beam shape provided rq and rmaxare independent of 80. In this case, 8, does not appear in the spd function and hence the value 60 = (a/2) can be used to simplify the problem. There is one case in which rq and rmaxand hence F(q) may depend on 60. This is when wall reflection effects are included. The effective beam width in the x-direction is then reduced by a factor depending on 60. In this case the beam can still be considered to be at a / 2 to the tube, provided the limits of integration are first calculated correctly as functions of do. Calculation of f ( q ) The method for finding f ( q ) described above is to determine F(q) using volume integrals and then differentiate with respect to q. Analytic integration can be difficult or sometimes impossible. The integrals, (2) and (3), often simplify in Cartesian coordinates. Analytic integration of only one integral is required since (2) is a special case of (3). It is sometimes advantageous to reverse the order of differentiation and integration by first differentiating under the integral in Vq thereby obtaining f (q) as f ( q ) = (dF/dq) = - ( d V q / d 4 ) ( l / V T ) .One is then led directly to the spd function without passing through F(q). This method is used in later sections of the paper. One disadvantage is that there are now two different integrals to evaluate. SPECTRAL POWERDENSITYFUNCTIONS Uniform Beam, Complete Insonation A . Rectangular Beam: The beam is assumed wider than the tube and it insonates a length L along the tube. The fraction of the insonated volume occupied by particles
__
ALDlS AND THOMPSON. CALCULATION OF DOPPLER DENSITY FUNCTIONS
1025
with velocities greater than q is
Therefore F(q) = 1 - (1 tion is f(q)
=
-
2- ( 1 n
q)2/"and hence the spd func-
- q)((2-")/")
(4)
This is the undistorted spectrum which could only be obtained under ideal conditions. The mean velocity is 4 = ( n / ( n + 2 ) ) , and 50% of the insonated volume is occupied by particles with velocities less than the median velocity qm = I - ( 1 / 2 > " / ~ . B. Circular Beam, Beam Radius A > a: The intersection of a circular beam and a circular tube is shown in Fig. 2(a). The total volume of intersection of the beam and the tube, obtained from ( 2 ) in Cartesian coordinates, is
t X
(a)
v'A?-x'
V, = 8 = 8
1:' j Sa
m
d
y
d
x
J ( A 2 - x 2 ) ( a 2- x2) dx.
Evaluation of this integral (Prudnikhov et al. [ 111) gives
t X
where h = a / A , and K and E are complete elliptic integrals of the first and second kinds, respectively, as defined in Abramowitz and Stegun [12]. The volume Vq is given immediately from (6) after replacing h by p = ( a / A ) (1 - 4)"". Therefore
F(q) = 1
-
(1 (1
+ P Z ) E ( P 2 ) - (1
+ h2)E(h*)
-
- P2)K(P2)
(1 - h2)K(h2).
(6)
(b)
Fig. 2 . (a) Region of intersection on the plane := 0 of a circular beam perpendicular to a circular tube (OD = 7r/2). The beam radius, A , is assumed larger than the tube radius a. (b) As in (a) but with wall reflection effects. The insonated width in the x-direction is reduced from 2a to 2 a ' .
In the limit as h 0 (7) reduces to the rectangular beam spd function, (4). +
The spd function is
Uniform Beam, Incomplete Insonation A. Rectangular Beam, Wall Rejection: Particles in a f(4) = n(1 - q ) ( ( l h 2 > ~ ( h-2 )(1 - h 2 ) ~ ( h 2 ) } ' region of flow which is not insonated are unable to con(7) tribute to the spd function. This incomplete insonation This was obtained after differentiating (6) by q and using could be due to a narrow rectangular beam or due to wall reflection. In Thompson et al. [14]it was shown that when the identities wall reflection effects are present, only that portion of the flow of width 2a ' in the x-direction (in the coordinates of Fig. 1 ) contributes to the spectrum where
+
3P2E(P2>
Equation (7) cannot be reduced to elementary functions, but the elliptic integrals are easily evaluated [ 131. When the tube and the beam are of equal radius h = 1 and (7) becomes
The angle 4,is the critical angle for the tube wall and the surrounding medium. Total external reflection of the beam occurs if the angle between the beam axis and the normal to the tube wall is greater than $ c . The spectrum shape does not depend on the Doppler angle OD in the case of a
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 39, NO. IO, OCTOBER 1992
1026
narrow rectangular beam, but it does when the incomplete insonation is due to wall reflection. Results are derived here for general a '. Once again we can perform the analysis with OD = 7r/2 since dependence on Doppler angle is represented by the value of a '. The total volume of the insonated region is
C. Circular Beam, Wall Rejection: With wall reflection only a width of tube 2 a ' in the x-direction can be insonated by any beam [see Fig. 2(b)]. If the beam radius A is less than or equal to a ' , then the spectrum is not affected by wall reflection and is given by (13). Spectra for the other two cases, a' < A Ia and A > a , can be found by application of the preceding results. Details for the A > a case are given below. It is more convenient to obtain f ( q ) by differentiating under the integral sign. The total volume VT is given by
where 1 = (a ' / a ) < 1 . Insonation is complete for velocities in the range 1 - I" to 1 , but incomplete for 0 5 q < 1 - 1". Therefore, F(q) is defined over a split domain and
where D(1) = sin-' 1 gives the spd function
+
I-.
{
(
)
1 (1 - q)l/" '
a'
J(a2 - x2>(A2 - x 2 ) k. 0
The volume Vq is evaluated according to whether or not the radius of the flow cylinder r = a ( l - q)l/" is greater than a '. We find
O I q < l - l "
B. Circular Beam, Beam Radius A < a: In the absence of wall reflection a circular beam which is narrower than a tube will give spectra which are independent of the Doppler angle. The total volume of intersection is calculated as before but with the roles of a and A reversed. Hence,
lo
Xmax
Vq = 8
VT=$((~
j
Differentiating this
sin-'
f(4)=
VT = 8
J ( r 2 - x 2 ) ( A 2 - x2) dx
(14)
where x,,, = a ' for 0 I q < 1 - 1" (incomplete insonation) and xmaX = r for 1 - I" I q I 1 (complete insonation) and 1 = a we find since F ( ~ =) 1 - vq/vT,
+j$)a!?($)-(l-$)K($))
(1 1) where h = a / A > 1 . The formula for F(q) is split into two regions of q by the value q = 1 - (1 / h " ) which corresponds to a flow cylinder with radius equal to A . From (3), (6) and (1 1) we find 1 E ( l - q)3/" * F(q) =
i'-
(( + $) ($)
C(h) 1
I--. ((1 h 3C(h)
l - (1
-
-$)
.
a1
(E$)),
+ P2)E(P2) - (1 - P2>K(P2)),
+
(15)
d(a2 - x2)(A2 - x2)
wherep = h(1 - q)l/" and C(h) = (1 [l/h2])E(1/h2) - (1 - [1/h2])K(1/h2). The spd function can be easily found by differentiation of (12) to be
0 I q < 1 - 1 ;
l - - -1 c q 5 l h" -
(12)
1027
ALDIS AND THOMPSON: CALCULATION OF DOPPLER DENSITY FUNCTIONS
Evaluation of the integrals [ 111 gives
0 5 q < l - Z "
-1
E(h2(1 - q ) 2 / " ) ,
where
~ ( hI ),
=
1 - 1"
Iq 5
1 (16)
+ h2)E(sin-' zlh2) - (1 - h 2 ) F(sin-' 1lh2) + lh2J(1 - z2h2)(1 - 1').
(1
F(q51m)and E(q5Jm)are incomplete elliptic integrals of the first and second kind respectively as defined in Wolfram ~31,
S,
-4
-2
0
2
4
Rad181 dlatance
Fig. 3 . Intensity against radial distance (normalized units) for axisymmetric Gaussian (solid line), jinc (dashed line) and sinc (dotted line) beams. The - 3 dB point of the central lobe for each curve occurs at beam radius F = a , the tube radius. Parameters are a / 2 o , = 0.588, 2Xa = 3.228 and Xa = I .389 for the Gaussian, jinc, and sinc beams, respectively.
6
F(6lm)
=
so
(1 - m sin2 e ) - ( I / z ) de,
6
E(q5lm) =
(1 - m sin2
The spd function is given here byf(q) = - (dW,/dq)/ W , which is equivalent to
de.
Nonuniform Beam Pro$les The spectral power distribution function F(q) can be calculated again using volume fractions but now all volumes are weighted by a beam profile function. In this section details of F(q) are suppressed and the spd functions are calculated by differentiating under the integral sign. Fig. 3 shows the beam intensity as a function of nonnalized radial distance for the axisymmetric Gaussian, jinc, and sinc beams. Let i denote the beam radius at which the -3 dB point of the central lobe occurs. The spread parameters in Fig. 3 have been chosen so that F = a , the tube radius. Apart from the weak sidelobes of the jinc and sinc beams the three beam profiles are very similar. All three beams extend radially to infinity, but acoustic intensity falls off rapidly with radial distance. If b(x, y) represents the beam profile in beam coordinates, then the weighted volume W4 of the flow cylinder of radius r occupied by particles with velocities greater than q = 1 - ( r / a ) "is given by
B(x) J
7 7 dr
where r = a(1 - q)"". The above formula (17) forf(q) could also have been used to calculate spectra for uniform beams, which have b(x, y) = 1. The finite size of the uniform beam would then need to be taken into account in the definition of B(x). Equation (15) is such an example, in which B(x) =
-. A . Gaussian Beam: We define an elliptical Gaussian beam by b(x, y) = exp (-[x2/2c7;] - Lv2/2a:]). The beamwidth in the x (cross-tube) direction is determined by the parameter 0'.This can be either the major or the minor axis of the elliptical beam. Adopting the approach outlined above, B(x) is found to be
Substituting B(x) into (17) the integral with respect to y cancels leaving = 8
~ ' B ( x ) ~ d r
where
The spectral power distribution function is F(q) = 1 W 4 / W , where W , is the weighted total insonated volume.
C
exp
(-$)
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. IO, OCTOBER 1992
1028
Evaluating the integrals [12, p. 3761 leads to
where the sine integral, Si ( x ) , is defined by Si (x) =
[' sin t dt. t
Jo
C. Sinc Beam: The third nonuniform beam pattern we consider is the axisymmetric sinc function
(18)
as found previously by Bascom and Cobbold [5]. Here Zo and ZIare modified Bessel functions. Hence, an elliptical Gaussian beam with either its major or minor axis aligned with the tube give a spectrum identical to that for a Gaussian beam with circular symmetry. Only the scale factor o1 which corresponds to the cross-stream direction enters into the final formula forf(q). B. Jinc Beam: In this section we consider the axisymmetric jinc beam profile
where again X is a parameter for the spread of the beam. The integrated beam function, B(x), is given by a
B(x) = 4X2x
1 P
2xr
Jo(t) dt.
0
The integral in (22) can be evaluated in terms of Struve functions but it is more efficient here to expand the Bessel function Jo in a power series. After substitution into (17) and reversing the order of integration the spd function is found to be
N
where the parameter A represents the spread of the beam. As X becomes smaller the main lobe of the beam becomes wider. The integrated beam profile B(x) is
B(x)
=
4
s
J:( A X2(x2
0
m
)
+ y2)
dY.
From Gradsteyn and Ryzhik [15, p. 7071 this can be evaluated as
B(x) =
Unfortunately, the integrals cannot be evaluated in terms of elementary functions.
2 HI(2hr) X3x2
~
where H l ( x ) is a Struve function as defined in [12]. The spd function can then be found from (17) which after a change of variable in each integral becomes P I
where CY = 2 X 4 l - q)'/" and = 2Xa. The two integrals can be evaluated in terms of the generalized hypergeometric series 2F3 [15, p. 7771. After simplification and using the series expansion for the sine integral [16, p. 531 we find the spd function is
RESULTSAND DISCUSSION Fig. 4 shows the spd functions obtained with complete insonation by a uniform rectangular beam ( 4 ) for axisymmetric velocity profiles of the type given by ( 1 ) . When n = 2 the flow is parabolic and the spd function is constant. As n increases the flow profile becomes blunter and the spd function approaches a spike at q = 1 as n --t 00 (plug flow). As n + 0 the spd function approaches a spike at q = 0. Since we have assumed throughout that the beam axis is centered in the tube, the spectral distortion due to incomplete insonation or beam nonuniformity always amounts to underrepresentation of the slower velocities near the wall. In the limit n + 00 there could be no spectral distortion. The differences between spd functions given by ( 4 ) and those calculated for other conditions are greater for smaller values of n. Since physiological flows are usually modeled with n I2, the following figures are all for n = 2. Equivalent curves for other n can be readily obtained from the general formulae derived in this paper. Fig. 5 shows spectra for a uniform circular beam with different ratios of the tube radius to the beam radius, h = a / A . The spectra are from (7), (8), and (13). As h 0 the beam radius can be regarded as approaching infinity and the spectrum is therefore the same as for a wide rectangular beam. Even when the beam is only twice as wide as the tube ( h = 0.5) the spectrum is very close to a constant. For higher values of h the regions of flow near the tube wall receive relatively low insonation and the spectra -+
f(4)=
1
CY
Si (CY) (CY)-
CY
~
CY
~
4 1 - 4)
sin
Si
(6)- 2
+ cos
CY
+ sin /3 + cos p ~
P
~
ALDIS A N D THOMPSON: CALCULATION OF DOPPLER DENSITY FUNCTIONS
1029
. 6 6
0.2-
0.0
0.2
0.6
0.4
0.8
0
1 .o
20
40
Doppkr mglo Bo (-1
q
80
80
(a)
Fig. 4 . Spectra for complete insonation by a uniform rectangular beam plotted against normalized velocity, q. The spectra are from (4) with different values of the parameter n in the velocity profile function (1). Parabolic flow (n = 2) is shown by the dashed line.
4
1
0 0.0
II I
I
0.2
0.4
I
I
0.6
0.8
1 .o
q
(b)
4 0.0
0.2
0.6
0.4
0.8
1 .o
q Fig. 5. Spectra for insonation by uniform circular beams of different radii with parabolic flow. The parameter h = a / A represents the ratio of tube radius to beam radius. The spectra are obtained from (7), (8), and (13).
show a greater weighting of high speed particles. Spectral distortion becomes a significant effect when the beam radius is less than the tube radius (h > 1). When the beam is narrower than the tube an infinite slope is observed at the point where the two sections of (1 3) meet. The effect of wall reflection on the Doppler spectrum has been discussed previously in Thompson et al. [ 141. It was shown that the effective beamwidth in one direction is reduced by an amount which depends on the Doppler angle and the critical angle [see (9)]. The relationship is shown in Fig. 6(a). It is apparent that distortion of the spectrum is more severe as BD decreases or as $, decreases. The effect for a uniform rectangular beam has been illustrated in [ 141. Fig. 6(b) shows spectra for a uniform circular beam with h = 1 and with wall reflection effects. A circular beam leads to some underrepresentation of low velocities, as does wall reflection. Wall reflection effects are thus less apparent with circular beams since the lower velocities are already underrepresented. This is particularly so when the beam is narrower than the tube, h > 1.
Fig. 6 . (a) The beamwidth reduction parameter I = a ’ / . plotted against Doppler angle Bo for several values of the critical angle IC., as given by (9). (b) Spectra for insonation by a uniform circular beam with h = 1 and with wall reflection present. Spectra for several values of the beamwidth reduction parameter 1 are shown. The parameter 1 combines the effects of Doppler angle and critical angle.
In the preceding section spd functions were derived for three nonuniform beam profiles: elliptical Gaussian, and axisymmetric jinc and sinc. The spectra from elliptical and axisymmetric Gaussian beams were shown to be equivalent provided one of the beam axes was aligned with the tube. The central lobes of the beams are very similar (see Fig. 3) when the beamwidth parameters are chosen as described to align the - 3 dB points. Consequently, the spd functions are very similar, except when the beam is narrow and the tube contains sidelobes. Fig. 7 shows the spd functions for the Gaussian, jinc and sinc beams with a / ? = 0.75, 1 and 2. In Fig. 8 the results for uniform circular beams (7), (8), and (13) and the Gaussian beam (18) are compared for different beamwidths. The radius, A , of the circular beam is set to the - 3 dB halfwidth of the Gaussian beam, i , so that h = a / A = a / ? . Spectra from the two beam types are again very similar provided the beams are wider than the tube. Spectra for a Gaussian beam are easily obtained from the single functional relationship, (1 8). The uniform circular beam is less physically realistic than the Gaussian beam. The spectra are not given by a single equation and show more distortion (see Fig. 8). Nonuniform beams are
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. IO, OCTOBER 1992
1030
0.24
0.0
0.2
0.4
0.6
0.8
1 .o
4 Fig. 7. Spectra for Gaussian (solid curves), jinc (dashed curves) and sin, (dotted curves) at a / ? = 0.75, 1, and 2 where i is the radius at which the -3 dB point occurs.
2.0
-*
I-
Gausslanl
a/P.2.0
,,7+
1.5
C
U al L
al
g 1.0 n
-2 e
0.5 U)
0.0 0.0
0.2
0.4
0.6
0.8
1
9
Fig. 8. Comparison of spectra from Gaussian (solid curves) and uniform circular beams (dashed curves) for various beamwidths. The circular beam radius A is set to the - 3 dB halfwidth of the Gaussian beam, i , so that h = a / A = a/?.
not much harder to work with analytically than uniform beams. Spectra for a jinc beam are particularly easy to evaluate since the formula, (20), involves only sines, cosines and the sine integral. The computations could be performed on a hand calculator. Spectra for the sinc beam, (23), involve integrals of Bessel functions and are more difficult to evaluate. The importance of these results is not so much in the numerical differences between the spectra studied (which are often small) but in the general applicability of the method given. For example, the effect of wall reflection on the spectra from nonuniform beams can be readily treated by changing limits of integration. Also the analysis which led to (2) and (3) was for a general beam shape A(0) and application to irregular beam shapes is possible. Even extreme or infolded beam shapes can be treated although the regions of integration could become quite complicated. Cases in which the beam does not pass
through the tube axis can also be treated. Formulae for the mean and median velocities can also be easily written down. A major assumption of the method described here is that the flow is axisymmetric. In fact some nonaxisymmetric cases, for example elliptical pipe flow, can also be treated. The flow was also assumed steady and of a power law form, (1). The steady flow assumption is not a practical restriction since the velocity profiles for many unsteady flows will only vary slightly with axial position over the region of insonation. Application to axisymmetric flow profiles other than (1) is straightforward. Provided the velocity is monotonic in r then VT and Vq are easily specified. When the velocity is not monotonic it will still be piecewise monotonic and spectra can be found using several regions of integration. In [14] and [lo] the spd functions were obtained using the chain rule of differentiation. That approach is only possible when the velocity is monotonic. Spectral distortion under the various conditions of incomplete and nonuniform insonation considered in this paper involves an overrepresentation of the higher velocity (axial) particles. Increasing the value of the parameter n in the velocity profile (1) leads to a similar change in spectral shape (Fig. 4). On the other hand, different flows may give rise to the same theoretical spd function. For example, the spd functions for plane shear flow and viscous flow through all tubes of elliptical cross section are identical. The observed Doppler spectrum may therefore be unchanged over a region of gradual elliptical compression of a circular tube. These underline the potential for error which exists if attempts are made to interpret velocity information from Doppler spectra without a good understanding of both the underlying flow and the insonating beam. REFERENCES [ l ] D. H. Evans, W. N. McDicken, R. Skidmore, and J . P. Woodcock, Doppler Ultrasound. Physics, Instrumentation, and Clinical Applications. New York: Wiley, 1989. [2] D. H. Evans, “Some aspects of the relationships between instantaneous volumetric blood flow and continuous wave Doppler ultrasound recordings-I. The effect of ultrasonic beam width on the output of maximum, mean and RMS frequency processors,” Ultrasound Med. Biol., vol. 8, pp. 605-609, 1982. [3] R. W. Gill, “Accuracy calculations for ultrasonic pulsed Doppler blood flow measurements,” Australas. Phys. Eng. Sci. Med., vol. 5 , pp. 51-57, 1982. [4] R. S . C. Cobbold, P. H. Veltink, and K. W. Johnston, “Influence of beam profile and degree of insonation on the CW Doppler ultrasound spectrum and mean velocity,” IEEE Trans. Son. Ultrason., vol. SU-30, pp. 364-370, 1983. [5] P. A. J. Bascom and R. S . C. Cobbold, “Effects of transducer beam geometry and flow velocity profile on the Doppler power spectrum: A theoretical study,” Ultrasound Med. Biol., vol. 16, pp. 279-295, 1990. [6] P. A. J. Bascom, R . S . C. Cobbold, and B. H. M .Roelofs, “Influence of spectral broadening on continuous wave Doppler ultrasound spectra: A geometric approach,” UltrasoundMed. Biol., vol. 12, pp. 387-395, 1986. [7] V . L. Newhouse, D. Censor, T. Vontz, J. A. Cisneros, and B. B. Goldberg, “Ultrasound Doppler probing of flows transverse with re-
ALDIS A N D THOMPSON: CALCULATION OF DOPPLER DENSITY FUNCTIONS
spect to the flow axis,’’ IEEE Trans. Biomed. Eng., vol. BME-34, pp. 779-789, 1987. [8] D. Censor, V. L. Newhouse, T. Vontz, and H. V. Ortega, “Theory of ultrasound Doppler-spectra velocimetry for arbitrary beam and flow configurations,” IEEE Trans. Biomed. Eng., vol. BME-35, pp. 740751, 1988. [9] G. L. Gooberman, Ultrasonics: Theory and Applications. London: The English Universities Press, 1968. [IO] J. M. J . G. Roevros, “Analogue processing of C. W. Doppler flowmeter signals to determine average frequency shift momentaneously without the use of a wave analyzer,” in Cardiovascular Applications of Ultrasound, R. S . Reneman and M. P. Spencer, Eds. Amsterdam: North-Holland, 1974, pp. 43-54. I l l ] A. P. Prudnikov, Y. A. Brychkov, and 0. I. Marichev, Integrals and Series. New York: Gordon and Breach, 1986. [I21 M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions. New York: Dover, 1972. [I31 S . Wolfram, Mathematica: A System for Doing Mathematics by Computer. Redwood City: Addison-Wesley, 1988. [I41 R. S . Thompson, G. K. Aldis, and I. W. Linnett, “Doppler ultrasound spectral power density distribution: Measurement artefacts in steady flow,” Med. B i d . Eng. Comput., vol. 28, pp. 60-66, 1990. [I51 I. S . Gradsteyn and I. W. Ryzhik, Table of Integrals, Series, and Products. New York: Academic, 1965, 4th ed. [I61 E. R. Hansen, A Table of Series and Products. Englewood Cliffs, NJ: Prentice-Hall, 1975.
1031
Geoffrey K. Aldis received the B.Sc. degree (Med.) from the University of Sydney in 1976, the B.Math. degree from the University of Wollongong in 1980, and the Ph.D. degree from the Department of Applied Mathematics and Theoretical Physics at the University of Cambridge in 1984 with work on the solute concentration fields near transporting epithelia. He is currently a lecturer in the Department of Mathematics at the University College, ADFA. His research interests include osmosis and physiological porous media flows. Rosemary S. Thompson received the B.Sc. (Hons 1) degree in applied mathematics and physics from the Australian National University, Canberra, in 1974; her interest in biomathematics led to the Ph.D. study of networks of model neurons in the Department of Applied Mathematics, University of Sydney (1982). She has worked on Doppler ultrasound techniques and their application to assessment of the fetal circulation. She is currently a Research Officer supported by NH&MRC grant with the University of Sydney at Westmead Hospital.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. IO. OCTOBER 1992
Calculation of Doppler Spectral Power Density Functions Geoffrey K. Aldis and Rosemary S . Thompson
Abstract-A volume integral method for the calculation of Doppler ultrasound spectral power density (spd) functions is described. Axisymmetric flow in a circular tube with a power law velocity profile is assumed. The spd function is regarded as a probability density function for scatterer velocity, and the assumptions under which this is justified are considered. It is shown that the spd function is independent of Doppler angle except in the presence of wall reflection effects. A coordinate system centered on the beam is used and this enables the integrals to be easily formulated for arbitrary beams. Irregularly shaped and nonuniform beams can be treated. For the common flow and beam patterns, which exhibit symmetry, the volume integrals can often be reduced to a single integral and evaluated directly. The method is applied and the spectra are calculated for various different cases. Results are obtained for uniform rectangular and circular insonating beams, and for nonuniform beams with Gaussian, jinc, and sinc profiles. The effects of narrow beams and wall reflection are shown. The method may be readily applied to other beam and flow patterns, and extension to more complicated situations is also discussed.
INTRODUCTION OPPLER spectra are readily obtained from different arteries in the body, but care is required in making inferences about physiological flow patterns from these spectra. The Doppler spectrum results from the flow velocity profile in the artery, but it may be modified by several other processes which are not flow-related. The spectral power density (the fraction of backscattered power contributed by scatterers at each velocity) may not be truly representative of the velocity density (the fraction of scatterers moving with each velocity). Following Evans et al. [l] the reasons for this can be summarized into three groups: i) scattering properties of blood or other target medium; ii) processing of the received signal by, for example, filtering or the FFT algorithm; iii) three-dimensional (3-D) geometry of the transduction process. In this paper the third factor, spectral distortion due to the 3-D geometry of the transduction process, is considered for continuous wave (CW) Doppler spectra. A num-
D
Manuscript received April 3, 1991; revised Februaty 10, 1992. This work was supported by the National Health and Medical Research Council of Australia. G. K. Aldis is with the Department of Mathematics, University College, University of New South Wales, Australian Defence Force Academy, Campbell ACT 2600, Australia. R. S . Thompson is with the Department of Obstetrics and Gynaecology, University of Sydney at Westmead Hospital, Westmead, NSW 2145, Australia. IEEE Log Number 9202575.
ber of different authors have previously discussed various aspects of this spectral distortion. Evans [2] and Evans et al. [ 11 considered the problem two-dimensionally using rectangular beams which may be narrower than the tube. Gill [3] examined errors in the estimated mean velocity, and presented some numerical calculations of spectra using a nonuniform beam profile. Cobbold et aZ. [4] developed a computer model of the 3-D problem incorporating effects due to a nonuniform beam profile, off-axis alignment of the beam, and differential attenuation between blood and tissue. Numerical results were obtained for various cases of interest and this approach was continued by Bascom and Cobbold [ 5 ] , who also derived an analytic expression for the spectral power density with a Gaussian beam profile. Cobbold’s model was also extended to study intrinsic spectral broadening by Bascom et al. [6]. Intrinsic spectral broadening is broadening of the received spectrum due to properties of the measurement system itself, and has been discussed extensively by other authors [l]. Newhouse et al. [7] and Censor et al. [8] presented theoretical results using the ultrasonic fields appropriate to certain transducers, whereas Bascom et al. [6] used a geometrical ray tracing approach. Spectral broadening can be described either geometrically or in terms of the finite transit time of scatterers through the beam [l]. Although it is a consequence of the 3-D geometry of the transduction process, spectral broadening is not considered here. This is equivalent to the assumption that all ultrasonic wavefronts are planar. Spectral broadening tends to smooth the shape of Doppler spectra [l] but since all frequencies are affected, this may not be evident except towards the extremities of the spectrum. The neglect of spectral broadening in this paper means that the calculated spectra will have an unphysical sharp cut off at the maximum frequency. In this paper we assume that the flow being insonated is steady axisymmetric flow of a simple power law form. We also assume that over the insonated region of the tube the beam has plane wavefronts and a fixed but not necessarily uniform cross-sectional intensity pattern. A circular transducer produces a beam for which the directivity function in the far field (Fraunhofer zone) is a jinc function, whereas for a finite length line source the directivity function is a sinc function [9]. The beam intensity pattern of ultrasonic transducers may be approximated using these
00 18-9294/92$03.OO 0 1992 IEEE
ALDIS AND THOMPSON: CALCULATION OF DOPPLER DENSITY FUNCTIONS
functions. A Gaussian beam intensity is also commonly used as it is simpler, and closely approximates both these functions in the axial region. Spectral power density (spd) functions can be calculated using volume fractions (Roevros [lo], Cobbold et al. [4]).A volume integral approach is developed in this paper using the 3-D geometry of the beam/tube configuration and an arbitrary beam profile. Integrals are formulated in beam coordinates which enables arbitrary beam profiles and shapes to be treated. It is shown that the spd function is independent of Doppler angle except in the presence of wall reflection effects. Even when the spd functions cannot be obtained analytically the volume integrals can often be simplified and computation reduced. A number of results are obtained using the method. Spectra are calculated for nonuniform Gaussian, jinc, and sinc beam profiles as well as for uniform rectangular and uniform circular beams. The effects of narrow beams and wall reflection are shown. METHODOF CALCULATION OF SPECTRAL POWER DENSITYFUNCTIONS
Assumptions The backscattered power from a small control volume of fluid is assumed proportional to the total number of scatterers in the volume and to the insonating beam intensity. We assume that there are a large number of independent point scatterers uniformly distributed throughout the insonated volume. Backscattered power is therefore proportional to volume if the insonating beam is uniform. When the beam has a variation in intensity over its crosssection the backscattered power is proportional to the volume weighted by the beam intensity. Over the insonated region of the tube the ultrasonic wavefronts are assumed perpendicular to the transducer face (no spectral broadening). There is no variation in the beam intensity pattern with axial distance from the transducer (thin-walled tube, no attenuation). We consider axisymmetric flow in a circular tube of radius a of the form
where q is the normalized velocity, r is radial distance from the tube axis, and n is a parameter. Poiseuille flow corresponds to n = 2 and as n -, 00 the profile approaches plug flow. All particles moving with velocity greater than q lie inside a flow cylinder of radius r = a ( l - q ) ' / " .The range of the normalized velocity q is 0 to 1 , but extension to any velocity range is easily accomplished by multiplying q by qmax.The Doppler spectra are presented in this study as functions of normalized velocity. The spectra could be equivalently regarded as functions of normalized Doppler frequency shift, since velocity and frequency shift are directly proportional by the Doppler equation.
1023
Spectral Power Distribution Function, F(q) Under the ideal conditions of uniform scatterer distribution and uniform insonation assumed here, the spd function is proportional to the velocity density and can be regarded as a probability density function for both spectral power and velocity. Under nonideal conditions the spd function may not be a probability density function for the velocity. Interest here is primarily on comparing the shape of observed spectra and not with the value of total backscattered power. The area under the spd function is therefore always left equal to 1 even under conditions of incomplete insonation. Associated with the spd function there i s a spectral probability distribution function, F(q) = j i f ( q ) dq. The value of the spectral power distribution function, F(q), is the fraction of the total insonated volume occupied by particles with velocities between 0 and q. F is a monotonic increasing function with F(0) = 0 and F( 1) = 1. It is often easier to determine F(q) and then to calculate the spd function from it by differentiation. The mean velocity, 4, can be obtained from Z j = Sh ( 1 - F(q)) dq = q f ( q ) dq. The median velocity qmis found from solving F(q,) = ( V I . The angle between the beam axis and the tube axis is the Doppler angle For plane wavefronts and axisym# metric flow a Doppler shift is only obtained for (a/2). However, it greatly simplifies the volume calculations if the transducer is inclined at a / 2 to the flow axis. We show below that for axisymmetric flow in a circular tube choosing Bo = (a/2) leads to the correct spd function. This holds even in the case of wall reflection, where the spd function depends on OD, provided the beam shape is appropriately modified. Calculation of F(q) Using Volume Integrals Integrals will be formulated using coordinate systems centered on the beam. The beam's cross-section has the general boundary A = A(%)and a uniform beam intensity over the beam cross-section is assumed for the moment. Fig. l(a) shows the beam and tube and their respective coordinate systems which share a common origin at the center of the tube. The tube Cartesian system (x', y', z ' ) is centered in the tube with the z '-axis along the tube axis and the y'-axis also in the plane of the page. The beam Cartesian system (x, y, z ) has its z-axis along the beam. The x- and x'-axes coincide. The y'- and z'-axes are obtained by rotating the y- and z-axes around the x-axis through OD. A cylindrical beam system ( r , 9, z) is also defined with x = r cos 9 and y = r sin 9. The surface of the tube is given by x f 2 y'* = a2 and since x ' = x and y ' = y cos - z sin 80 we find that the intersection of the beam and the tube occurs at
+
Fig. l(b) shows the projection of the tube onto the plane = 0 (a cross-section of the beam). The projection shows
z
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 39, NO. IO, OCTOBER 1992
1024
Hence
ir i
2
VT = sin 80 2'
rmax
r(a2 - r2 cos2
dr d€J ( 2 )
in which the dependence on Doppler angle is conveniently outside the integrals. Particles with velocities greater than q lie inside the flow cylinder of radius r = a ( l - q)'/".The volume occupied by the intersection of the beam with these particles can be obtained readily from the expression above for VTas
jr 1 r(a2(1
2 Vq = sin 60
'4
-
q)2/"- r2 cos2
dr de
(3) where
\
x=x'
X
(b) Fig. 1. (a) Cartesian coordinate systems for the beam ( x , y. z) and the tube (x', y', z ' ) . The tube system is obtained by rotating the beam system through Bo around the common x = 1'-axis. (b) A cross-section of the beam (in the plane z = 0) with the projection o f the tube shown shaded. Distance from the origin to the beam edge is A ( @ , distance to the projected tube edge is a / ( c o s 01.
the region of r and 0 which is used to calculate the total volume of intersection of the beam and tube. At each value of 0 integration extends in the r direction to either the edge of the tube, a l l c o s 81, or to the edge of the beam, A ( @ , whichever is encountered first. The total volume of the intersection of the beam and the tube, V,, is therefore
where
- (a2 -
r2 cos2
+ ( a 2 - r2 cos2
The volume fraction occupied by particles with velocities greater than q is V q / V p Therefore F(4), the volume fraction occupied by particles with velocities less than q , is F(q) = 1 - V q / V T . The dependence on Doppler angle cancels in F(q) for any beam shape provided rq and rmaxare independent of 80. In this case, 8, does not appear in the spd function and hence the value 60 = (a/2) can be used to simplify the problem. There is one case in which rq and rmaxand hence F(q) may depend on 60. This is when wall reflection effects are included. The effective beam width in the x-direction is then reduced by a factor depending on 60. In this case the beam can still be considered to be at a / 2 to the tube, provided the limits of integration are first calculated correctly as functions of do. Calculation of f ( q ) The method for finding f ( q ) described above is to determine F(q) using volume integrals and then differentiate with respect to q. Analytic integration can be difficult or sometimes impossible. The integrals, (2) and (3), often simplify in Cartesian coordinates. Analytic integration of only one integral is required since (2) is a special case of (3). It is sometimes advantageous to reverse the order of differentiation and integration by first differentiating under the integral in Vq thereby obtaining f (q) as f ( q ) = (dF/dq) = - ( d V q / d 4 ) ( l / V T ) .One is then led directly to the spd function without passing through F(q). This method is used in later sections of the paper. One disadvantage is that there are now two different integrals to evaluate. SPECTRAL POWERDENSITYFUNCTIONS Uniform Beam, Complete Insonation A . Rectangular Beam: The beam is assumed wider than the tube and it insonates a length L along the tube. The fraction of the insonated volume occupied by particles
__
ALDlS AND THOMPSON. CALCULATION OF DOPPLER DENSITY FUNCTIONS
1025
with velocities greater than q is
Therefore F(q) = 1 - (1 tion is f(q)
=
-
2- ( 1 n
q)2/"and hence the spd func-
- q)((2-")/")
(4)
This is the undistorted spectrum which could only be obtained under ideal conditions. The mean velocity is 4 = ( n / ( n + 2 ) ) , and 50% of the insonated volume is occupied by particles with velocities less than the median velocity qm = I - ( 1 / 2 > " / ~ . B. Circular Beam, Beam Radius A > a: The intersection of a circular beam and a circular tube is shown in Fig. 2(a). The total volume of intersection of the beam and the tube, obtained from ( 2 ) in Cartesian coordinates, is
t X
(a)
v'A?-x'
V, = 8 = 8
1:' j Sa
m
d
y
d
x
J ( A 2 - x 2 ) ( a 2- x2) dx.
Evaluation of this integral (Prudnikhov et al. [ 111) gives
t X
where h = a / A , and K and E are complete elliptic integrals of the first and second kinds, respectively, as defined in Abramowitz and Stegun [12]. The volume Vq is given immediately from (6) after replacing h by p = ( a / A ) (1 - 4)"". Therefore
F(q) = 1
-
(1 (1
+ P Z ) E ( P 2 ) - (1
+ h2)E(h*)
-
- P2)K(P2)
(1 - h2)K(h2).
(6)
(b)
Fig. 2 . (a) Region of intersection on the plane := 0 of a circular beam perpendicular to a circular tube (OD = 7r/2). The beam radius, A , is assumed larger than the tube radius a. (b) As in (a) but with wall reflection effects. The insonated width in the x-direction is reduced from 2a to 2 a ' .
In the limit as h 0 (7) reduces to the rectangular beam spd function, (4). +
The spd function is
Uniform Beam, Incomplete Insonation A. Rectangular Beam, Wall Rejection: Particles in a f(4) = n(1 - q ) ( ( l h 2 > ~ ( h-2 )(1 - h 2 ) ~ ( h 2 ) } ' region of flow which is not insonated are unable to con(7) tribute to the spd function. This incomplete insonation This was obtained after differentiating (6) by q and using could be due to a narrow rectangular beam or due to wall reflection. In Thompson et al. [14]it was shown that when the identities wall reflection effects are present, only that portion of the flow of width 2a ' in the x-direction (in the coordinates of Fig. 1 ) contributes to the spectrum where
+
3P2E(P2>
Equation (7) cannot be reduced to elementary functions, but the elliptic integrals are easily evaluated [ 131. When the tube and the beam are of equal radius h = 1 and (7) becomes
The angle 4,is the critical angle for the tube wall and the surrounding medium. Total external reflection of the beam occurs if the angle between the beam axis and the normal to the tube wall is greater than $ c . The spectrum shape does not depend on the Doppler angle OD in the case of a
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 39, NO. IO, OCTOBER 1992
1026
narrow rectangular beam, but it does when the incomplete insonation is due to wall reflection. Results are derived here for general a '. Once again we can perform the analysis with OD = 7r/2 since dependence on Doppler angle is represented by the value of a '. The total volume of the insonated region is
C. Circular Beam, Wall Rejection: With wall reflection only a width of tube 2 a ' in the x-direction can be insonated by any beam [see Fig. 2(b)]. If the beam radius A is less than or equal to a ' , then the spectrum is not affected by wall reflection and is given by (13). Spectra for the other two cases, a' < A Ia and A > a , can be found by application of the preceding results. Details for the A > a case are given below. It is more convenient to obtain f ( q ) by differentiating under the integral sign. The total volume VT is given by
where 1 = (a ' / a ) < 1 . Insonation is complete for velocities in the range 1 - I" to 1 , but incomplete for 0 5 q < 1 - 1". Therefore, F(q) is defined over a split domain and
where D(1) = sin-' 1 gives the spd function
+
I-.
{
(
)
1 (1 - q)l/" '
a'
J(a2 - x2>(A2 - x 2 ) k. 0
The volume Vq is evaluated according to whether or not the radius of the flow cylinder r = a ( l - q)l/" is greater than a '. We find
O I q < l - l "
B. Circular Beam, Beam Radius A < a: In the absence of wall reflection a circular beam which is narrower than a tube will give spectra which are independent of the Doppler angle. The total volume of intersection is calculated as before but with the roles of a and A reversed. Hence,
lo
Xmax
Vq = 8
VT=$((~
j
Differentiating this
sin-'
f(4)=
VT = 8
J ( r 2 - x 2 ) ( A 2 - x2) dx
(14)
where x,,, = a ' for 0 I q < 1 - 1" (incomplete insonation) and xmaX = r for 1 - I" I q I 1 (complete insonation) and 1 = a we find since F ( ~ =) 1 - vq/vT,
+j$)a!?($)-(l-$)K($))
(1 1) where h = a / A > 1 . The formula for F(q) is split into two regions of q by the value q = 1 - (1 / h " ) which corresponds to a flow cylinder with radius equal to A . From (3), (6) and (1 1) we find 1 E ( l - q)3/" * F(q) =
i'-
(( + $) ($)
C(h) 1
I--. ((1 h 3C(h)
l - (1
-
-$)
.
a1
(E$)),
+ P2)E(P2) - (1 - P2>K(P2)),
+
(15)
d(a2 - x2)(A2 - x2)
wherep = h(1 - q)l/" and C(h) = (1 [l/h2])E(1/h2) - (1 - [1/h2])K(1/h2). The spd function can be easily found by differentiation of (12) to be
0 I q < 1 - 1 ;
l - - -1 c q 5 l h" -
(12)
1027
ALDIS AND THOMPSON: CALCULATION OF DOPPLER DENSITY FUNCTIONS
Evaluation of the integrals [ 111 gives
0 5 q < l - Z "
-1
E(h2(1 - q ) 2 / " ) ,
where
~ ( hI ),
=
1 - 1"
Iq 5
1 (16)
+ h2)E(sin-' zlh2) - (1 - h 2 ) F(sin-' 1lh2) + lh2J(1 - z2h2)(1 - 1').
(1
F(q51m)and E(q5Jm)are incomplete elliptic integrals of the first and second kind respectively as defined in Wolfram ~31,
S,
-4
-2
0
2
4
Rad181 dlatance
Fig. 3 . Intensity against radial distance (normalized units) for axisymmetric Gaussian (solid line), jinc (dashed line) and sinc (dotted line) beams. The - 3 dB point of the central lobe for each curve occurs at beam radius F = a , the tube radius. Parameters are a / 2 o , = 0.588, 2Xa = 3.228 and Xa = I .389 for the Gaussian, jinc, and sinc beams, respectively.
6
F(6lm)
=
so
(1 - m sin2 e ) - ( I / z ) de,
6
E(q5lm) =
(1 - m sin2
The spd function is given here byf(q) = - (dW,/dq)/ W , which is equivalent to
de.
Nonuniform Beam Pro$les The spectral power distribution function F(q) can be calculated again using volume fractions but now all volumes are weighted by a beam profile function. In this section details of F(q) are suppressed and the spd functions are calculated by differentiating under the integral sign. Fig. 3 shows the beam intensity as a function of nonnalized radial distance for the axisymmetric Gaussian, jinc, and sinc beams. Let i denote the beam radius at which the -3 dB point of the central lobe occurs. The spread parameters in Fig. 3 have been chosen so that F = a , the tube radius. Apart from the weak sidelobes of the jinc and sinc beams the three beam profiles are very similar. All three beams extend radially to infinity, but acoustic intensity falls off rapidly with radial distance. If b(x, y) represents the beam profile in beam coordinates, then the weighted volume W4 of the flow cylinder of radius r occupied by particles with velocities greater than q = 1 - ( r / a ) "is given by
B(x) J
7 7 dr
where r = a(1 - q)"". The above formula (17) forf(q) could also have been used to calculate spectra for uniform beams, which have b(x, y) = 1. The finite size of the uniform beam would then need to be taken into account in the definition of B(x). Equation (15) is such an example, in which B(x) =
-. A . Gaussian Beam: We define an elliptical Gaussian beam by b(x, y) = exp (-[x2/2c7;] - Lv2/2a:]). The beamwidth in the x (cross-tube) direction is determined by the parameter 0'.This can be either the major or the minor axis of the elliptical beam. Adopting the approach outlined above, B(x) is found to be
Substituting B(x) into (17) the integral with respect to y cancels leaving = 8
~ ' B ( x ) ~ d r
where
The spectral power distribution function is F(q) = 1 W 4 / W , where W , is the weighted total insonated volume.
C
exp
(-$)
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. IO, OCTOBER 1992
1028
Evaluating the integrals [12, p. 3761 leads to
where the sine integral, Si ( x ) , is defined by Si (x) =
[' sin t dt. t
Jo
C. Sinc Beam: The third nonuniform beam pattern we consider is the axisymmetric sinc function
(18)
as found previously by Bascom and Cobbold [5]. Here Zo and ZIare modified Bessel functions. Hence, an elliptical Gaussian beam with either its major or minor axis aligned with the tube give a spectrum identical to that for a Gaussian beam with circular symmetry. Only the scale factor o1 which corresponds to the cross-stream direction enters into the final formula forf(q). B. Jinc Beam: In this section we consider the axisymmetric jinc beam profile
where again X is a parameter for the spread of the beam. The integrated beam function, B(x), is given by a
B(x) = 4X2x
1 P
2xr
Jo(t) dt.
0
The integral in (22) can be evaluated in terms of Struve functions but it is more efficient here to expand the Bessel function Jo in a power series. After substitution into (17) and reversing the order of integration the spd function is found to be
N
where the parameter A represents the spread of the beam. As X becomes smaller the main lobe of the beam becomes wider. The integrated beam profile B(x) is
B(x)
=
4
s
J:( A X2(x2
0
m
)
+ y2)
dY.
From Gradsteyn and Ryzhik [15, p. 7071 this can be evaluated as
B(x) =
Unfortunately, the integrals cannot be evaluated in terms of elementary functions.
2 HI(2hr) X3x2
~
where H l ( x ) is a Struve function as defined in [12]. The spd function can then be found from (17) which after a change of variable in each integral becomes P I
where CY = 2 X 4 l - q)'/" and = 2Xa. The two integrals can be evaluated in terms of the generalized hypergeometric series 2F3 [15, p. 7771. After simplification and using the series expansion for the sine integral [16, p. 531 we find the spd function is
RESULTSAND DISCUSSION Fig. 4 shows the spd functions obtained with complete insonation by a uniform rectangular beam ( 4 ) for axisymmetric velocity profiles of the type given by ( 1 ) . When n = 2 the flow is parabolic and the spd function is constant. As n increases the flow profile becomes blunter and the spd function approaches a spike at q = 1 as n --t 00 (plug flow). As n + 0 the spd function approaches a spike at q = 0. Since we have assumed throughout that the beam axis is centered in the tube, the spectral distortion due to incomplete insonation or beam nonuniformity always amounts to underrepresentation of the slower velocities near the wall. In the limit n + 00 there could be no spectral distortion. The differences between spd functions given by ( 4 ) and those calculated for other conditions are greater for smaller values of n. Since physiological flows are usually modeled with n I2, the following figures are all for n = 2. Equivalent curves for other n can be readily obtained from the general formulae derived in this paper. Fig. 5 shows spectra for a uniform circular beam with different ratios of the tube radius to the beam radius, h = a / A . The spectra are from (7), (8), and (13). As h 0 the beam radius can be regarded as approaching infinity and the spectrum is therefore the same as for a wide rectangular beam. Even when the beam is only twice as wide as the tube ( h = 0.5) the spectrum is very close to a constant. For higher values of h the regions of flow near the tube wall receive relatively low insonation and the spectra -+
f(4)=
1
CY
Si (CY) (CY)-
CY
~
CY
~
4 1 - 4)
sin
Si
(6)- 2
+ cos
CY
+ sin /3 + cos p ~
P
~
ALDIS A N D THOMPSON: CALCULATION OF DOPPLER DENSITY FUNCTIONS
1029
. 6 6
0.2-
0.0
0.2
0.6
0.4
0.8
0
1 .o
20
40
Doppkr mglo Bo (-1
q
80
80
(a)
Fig. 4 . Spectra for complete insonation by a uniform rectangular beam plotted against normalized velocity, q. The spectra are from (4) with different values of the parameter n in the velocity profile function (1). Parabolic flow (n = 2) is shown by the dashed line.
4
1
0 0.0
II I
I
0.2
0.4
I
I
0.6
0.8
1 .o
q
(b)
4 0.0
0.2
0.6
0.4
0.8
1 .o
q Fig. 5. Spectra for insonation by uniform circular beams of different radii with parabolic flow. The parameter h = a / A represents the ratio of tube radius to beam radius. The spectra are obtained from (7), (8), and (13).
show a greater weighting of high speed particles. Spectral distortion becomes a significant effect when the beam radius is less than the tube radius (h > 1). When the beam is narrower than the tube an infinite slope is observed at the point where the two sections of (1 3) meet. The effect of wall reflection on the Doppler spectrum has been discussed previously in Thompson et al. [ 141. It was shown that the effective beamwidth in one direction is reduced by an amount which depends on the Doppler angle and the critical angle [see (9)]. The relationship is shown in Fig. 6(a). It is apparent that distortion of the spectrum is more severe as BD decreases or as $, decreases. The effect for a uniform rectangular beam has been illustrated in [ 141. Fig. 6(b) shows spectra for a uniform circular beam with h = 1 and with wall reflection effects. A circular beam leads to some underrepresentation of low velocities, as does wall reflection. Wall reflection effects are thus less apparent with circular beams since the lower velocities are already underrepresented. This is particularly so when the beam is narrower than the tube, h > 1.
Fig. 6 . (a) The beamwidth reduction parameter I = a ’ / . plotted against Doppler angle Bo for several values of the critical angle IC., as given by (9). (b) Spectra for insonation by a uniform circular beam with h = 1 and with wall reflection present. Spectra for several values of the beamwidth reduction parameter 1 are shown. The parameter 1 combines the effects of Doppler angle and critical angle.
In the preceding section spd functions were derived for three nonuniform beam profiles: elliptical Gaussian, and axisymmetric jinc and sinc. The spectra from elliptical and axisymmetric Gaussian beams were shown to be equivalent provided one of the beam axes was aligned with the tube. The central lobes of the beams are very similar (see Fig. 3) when the beamwidth parameters are chosen as described to align the - 3 dB points. Consequently, the spd functions are very similar, except when the beam is narrow and the tube contains sidelobes. Fig. 7 shows the spd functions for the Gaussian, jinc and sinc beams with a / ? = 0.75, 1 and 2. In Fig. 8 the results for uniform circular beams (7), (8), and (13) and the Gaussian beam (18) are compared for different beamwidths. The radius, A , of the circular beam is set to the - 3 dB halfwidth of the Gaussian beam, i , so that h = a / A = a / ? . Spectra from the two beam types are again very similar provided the beams are wider than the tube. Spectra for a Gaussian beam are easily obtained from the single functional relationship, (1 8). The uniform circular beam is less physically realistic than the Gaussian beam. The spectra are not given by a single equation and show more distortion (see Fig. 8). Nonuniform beams are
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. IO, OCTOBER 1992
1030
0.24
0.0
0.2
0.4
0.6
0.8
1 .o
4 Fig. 7. Spectra for Gaussian (solid curves), jinc (dashed curves) and sin, (dotted curves) at a / ? = 0.75, 1, and 2 where i is the radius at which the -3 dB point occurs.
2.0
-*
I-
Gausslanl
a/P.2.0
,,7+
1.5
C
U al L
al
g 1.0 n
-2 e
0.5 U)
0.0 0.0
0.2
0.4
0.6
0.8
1
9
Fig. 8. Comparison of spectra from Gaussian (solid curves) and uniform circular beams (dashed curves) for various beamwidths. The circular beam radius A is set to the - 3 dB halfwidth of the Gaussian beam, i , so that h = a / A = a/?.
not much harder to work with analytically than uniform beams. Spectra for a jinc beam are particularly easy to evaluate since the formula, (20), involves only sines, cosines and the sine integral. The computations could be performed on a hand calculator. Spectra for the sinc beam, (23), involve integrals of Bessel functions and are more difficult to evaluate. The importance of these results is not so much in the numerical differences between the spectra studied (which are often small) but in the general applicability of the method given. For example, the effect of wall reflection on the spectra from nonuniform beams can be readily treated by changing limits of integration. Also the analysis which led to (2) and (3) was for a general beam shape A(0) and application to irregular beam shapes is possible. Even extreme or infolded beam shapes can be treated although the regions of integration could become quite complicated. Cases in which the beam does not pass
through the tube axis can also be treated. Formulae for the mean and median velocities can also be easily written down. A major assumption of the method described here is that the flow is axisymmetric. In fact some nonaxisymmetric cases, for example elliptical pipe flow, can also be treated. The flow was also assumed steady and of a power law form, (1). The steady flow assumption is not a practical restriction since the velocity profiles for many unsteady flows will only vary slightly with axial position over the region of insonation. Application to axisymmetric flow profiles other than (1) is straightforward. Provided the velocity is monotonic in r then VT and Vq are easily specified. When the velocity is not monotonic it will still be piecewise monotonic and spectra can be found using several regions of integration. In [14] and [lo] the spd functions were obtained using the chain rule of differentiation. That approach is only possible when the velocity is monotonic. Spectral distortion under the various conditions of incomplete and nonuniform insonation considered in this paper involves an overrepresentation of the higher velocity (axial) particles. Increasing the value of the parameter n in the velocity profile (1) leads to a similar change in spectral shape (Fig. 4). On the other hand, different flows may give rise to the same theoretical spd function. For example, the spd functions for plane shear flow and viscous flow through all tubes of elliptical cross section are identical. The observed Doppler spectrum may therefore be unchanged over a region of gradual elliptical compression of a circular tube. These underline the potential for error which exists if attempts are made to interpret velocity information from Doppler spectra without a good understanding of both the underlying flow and the insonating beam. REFERENCES [ l ] D. H. Evans, W. N. McDicken, R. Skidmore, and J . P. Woodcock, Doppler Ultrasound. Physics, Instrumentation, and Clinical Applications. New York: Wiley, 1989. [2] D. H. Evans, “Some aspects of the relationships between instantaneous volumetric blood flow and continuous wave Doppler ultrasound recordings-I. The effect of ultrasonic beam width on the output of maximum, mean and RMS frequency processors,” Ultrasound Med. Biol., vol. 8, pp. 605-609, 1982. [3] R. W. Gill, “Accuracy calculations for ultrasonic pulsed Doppler blood flow measurements,” Australas. Phys. Eng. Sci. Med., vol. 5 , pp. 51-57, 1982. [4] R. S . C. Cobbold, P. H. Veltink, and K. W. Johnston, “Influence of beam profile and degree of insonation on the CW Doppler ultrasound spectrum and mean velocity,” IEEE Trans. Son. Ultrason., vol. SU-30, pp. 364-370, 1983. [5] P. A. J. Bascom and R. S . C. Cobbold, “Effects of transducer beam geometry and flow velocity profile on the Doppler power spectrum: A theoretical study,” Ultrasound Med. Biol., vol. 16, pp. 279-295, 1990. [6] P. A. J. Bascom, R . S . C. Cobbold, and B. H. M .Roelofs, “Influence of spectral broadening on continuous wave Doppler ultrasound spectra: A geometric approach,” UltrasoundMed. Biol., vol. 12, pp. 387-395, 1986. [7] V . L. Newhouse, D. Censor, T. Vontz, J. A. Cisneros, and B. B. Goldberg, “Ultrasound Doppler probing of flows transverse with re-
ALDIS A N D THOMPSON: CALCULATION OF DOPPLER DENSITY FUNCTIONS
spect to the flow axis,’’ IEEE Trans. Biomed. Eng., vol. BME-34, pp. 779-789, 1987. [8] D. Censor, V. L. Newhouse, T. Vontz, and H. V. Ortega, “Theory of ultrasound Doppler-spectra velocimetry for arbitrary beam and flow configurations,” IEEE Trans. Biomed. Eng., vol. BME-35, pp. 740751, 1988. [9] G. L. Gooberman, Ultrasonics: Theory and Applications. London: The English Universities Press, 1968. [IO] J. M. J . G. Roevros, “Analogue processing of C. W. Doppler flowmeter signals to determine average frequency shift momentaneously without the use of a wave analyzer,” in Cardiovascular Applications of Ultrasound, R. S . Reneman and M. P. Spencer, Eds. Amsterdam: North-Holland, 1974, pp. 43-54. I l l ] A. P. Prudnikov, Y. A. Brychkov, and 0. I. Marichev, Integrals and Series. New York: Gordon and Breach, 1986. [I21 M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions. New York: Dover, 1972. [I31 S . Wolfram, Mathematica: A System for Doing Mathematics by Computer. Redwood City: Addison-Wesley, 1988. [I41 R. S . Thompson, G. K. Aldis, and I. W. Linnett, “Doppler ultrasound spectral power density distribution: Measurement artefacts in steady flow,” Med. B i d . Eng. Comput., vol. 28, pp. 60-66, 1990. [I51 I. S . Gradsteyn and I. W. Ryzhik, Table of Integrals, Series, and Products. New York: Academic, 1965, 4th ed. [I61 E. R. Hansen, A Table of Series and Products. Englewood Cliffs, NJ: Prentice-Hall, 1975.
1031
Geoffrey K. Aldis received the B.Sc. degree (Med.) from the University of Sydney in 1976, the B.Math. degree from the University of Wollongong in 1980, and the Ph.D. degree from the Department of Applied Mathematics and Theoretical Physics at the University of Cambridge in 1984 with work on the solute concentration fields near transporting epithelia. He is currently a lecturer in the Department of Mathematics at the University College, ADFA. His research interests include osmosis and physiological porous media flows. Rosemary S. Thompson received the B.Sc. (Hons 1) degree in applied mathematics and physics from the Australian National University, Canberra, in 1974; her interest in biomathematics led to the Ph.D. study of networks of model neurons in the Department of Applied Mathematics, University of Sydney (1982). She has worked on Doppler ultrasound techniques and their application to assessment of the fetal circulation. She is currently a Research Officer supported by NH&MRC grant with the University of Sydney at Westmead Hospital.
