I E E E TRANSACTIONS ON B I O M E D I C A L E N G I N E E R I N G , VOL. 38. NO. 6. J U N E 1991
589
Modern Spectral Analysis Techniques for Blood Flow Velocity and Spectral Measurements with Pulsed Doppler Ultrasound Jean-Yves David, Steven A. Jones, Member, IEEE, and Don P. Giddens
Abstract-Four spectral analysis techniques were applied to pulsed Doppler ultrasonic quadrature signals to compare the relative merits of each technique for estimation of flow velocity and Doppler spectra. The four techniques were 1) the fast Fourier transform method, 2) the maximum likelihood method, 3) the Burg autoregressive algorithm, and 4) the modified covariance approach to autoregressive modeling. Both simulated signals and signals obtained from an in vitro flow system were studied. Optimal parameter values (e.g., model orders) were determined for each method, and the effects of signal-to-noise ratio and signal bandwidth were investigated. The modern spectral analysis techniques were shown to be superior to Fourier techniques in most circumstances, provided the model order was chosen appropriately. Robustness considerations tended to recommend the maximum likelihood method for both velocity and spectral estimation. Despite the restrictions of steady laminar flow,, the results provide important basic information concerning the applicability of modern spectral analysis techniques to Doppler ultrasonic evaluation of arterial disease.
I . INTRODUCTION
P
ULSED Doppler ultrasound instruments have been used in clinical applications which range from the study of fetal blood flow to the diagnosis of atherosclerosis. In one of the most common signal analysis methods used with these instruments, the received pulses are downmixed to quadrature signals and these are transformed into power spectra. Although the principle of Doppler measurement is simple, the spectral content of quadrature signals is complicated by a number of devicerelated factors which introduce uncertainty in the velocity estimations. Low spatial resolution in regions of large velocity gradients, low temporal resolution for rapidly changing flows, instrument noise, and transit time effects all contribute to a broadening of the Doppler power spectra [ I ] . However, such broadening may also arise from limitations in the resolution and accuracy of the spectrum estimator employed. It is thus important to determine the appropriateness of the available spectrum estimator models to Doppler ultrasound signals. Manuscript received December 7. 1989: rebised August 13, 1990. J:Y. David is with the School o f Mechanical Engineering, Georgia Institute o f Technology. Atlanta. C A 30332. S . A. Jones and D. P. Giddens are with thc School o l Aerospace Engineering, Georgia Institute of Technology, Atlanta. C A 30332. IEEE Log Number 9144684.
A number of spectral estimation techniques have recently been developed and have been compared to the more standard fast Fourier transform (FFT) method for Doppler ultrasonic signal processing. The most common of these methods are autoregressive (AR), moving average, autoregressive moving average, and the maximum likelihood method (MLM). In addition, AR spectra can be computed via several algorithms. These include the Burg method, the Yule-Walker method, and the modified covariance method (MCM). These techniques are described in more detail by Kay and Marple [2]. D’luna et al. [3], [4] used the Burg method for AR modeling to investigate in virro poststenotic vortices. The final prediction error criterion (FPE) of Akaike [ 5 ] was used to determine the appropriate model order. Use of the AR was motivated by its greater spectral resolution, but no comparison of the autoregressive approach with other methods was attempted. Kaluzynski [6] compared FFT and AR methods and showed that for short data records, the AR yields more statistical stability than the FFT. He did not address the accuracy of the velocity estimation, but considered only the stability of the spectral shape. Again, the FPE was used for the model order selection. Talhami and Kitney [7] investigated an autoregressive moving average model but did not address the problem of model order selection. Moreover, they used simulated signals that did not show the random phase fluctuations characteristic of real Doppler signals. An extensive study of several modem spectral analysis methods, including AR, MLM, and autoregressive moving average, was done by Vaitkus and Cobbold [8], [9]. Simulated signals were used which modeled the shape of a spectrum taken in vivo with continuous wave ultrasound, but the simulation was not based on the physics of the Doppler process. The study was directed toward spectral estimation and did not discuss velocity estimation. Limited in vivo data were presented, and although the actual in vivo flow field was unknown, many of the results of the simulation studies were supported by these data. For velocity estimation, Jones and Giddens [ 101 provided a comparison between FFT methods and Burg AR modeling. This work can be viewed as an extension of these investigations. Few of the investigations to date have quantitatively examined the accuracy of the modern spectral analysis
0018-929419 110600-0589$01.OO O 1991 IEEE
590
I E E E T R A N S A C T I O N S O N BlOMEDlCAL E N G I N E E R I N G , VOL. 38. N O . 6, J U N E 1991
techniques. Those investigations which were quantitative have been undertaken with simulations rather than actual flow data. These simulations were designed to model spectral shape or amplitude modulation phenomena rather than the underlying physics of Doppler signals. Although much work has been done to examine spectral estimation, velocity estimation methods with modern techniques have not been extensively studied. The effect of signal-to-noise ratio on the accuracy of velocity measurements has also been neglected. This paper compares the performance of the FFT, MLM, MCM, and Burg AR under a variety of signall noise and signal bandwidth conditions. Simulated signals which incorporate the essential physics of Doppler measurement and signals derived from laboratory flow studies will be used for this purpose. Both spectral and velocity estimates will be examined. The spectral analysis methods to be studied were selected on the following basis. The Burg AR technique is known to perform well with signals of short duration; MCM exhibits globally better behavior than Burg AR for noisy sinusoids; MLM has a smaller variance than the autoregressive methods under some conditions [ 1 11, [ 121. The FFT was used for comparison since it is the technique commonly employed in commercially available devices. 11. METHODS A . Spectral Analysis Methods The spectral analysis methods employed have been described extensively in the literature (21, [121. The periodogram was used with either a rectangular or a Hanning window. Autoregressive spectra were obtained via both the Burg algorithm and the MCM approach. MLM spectra were evaluated by the method of Musicus [ 131, which uses the AR parameters from the Burg algorithm to estimate the spectrum. Each algorithm is implemented for complex signals in order to retain the directionality information of the Doppler quadrature signals. The convergence criteria for the MCM were switched off during the following comparisons so that the program could iterate until a specified order had been reached or until it terminated automatically upon detection of numerical difficulty. A detailed description of each estimator may be found in 121 and 112).
modulation is manifested by random fluctuations in the detected Doppler shifted frequency that are not related to actual flow fluctuations. This phenomenon has been termed “Doppler ambiguity” [ l ] , and it affects the velocity and spectral estimations. Briefly, the simulation model generates two quadrature signals expressed as
Q l ( t ) = A(t) cos
(U(,[)
+ B ( f )sin (w,,t)
(1)
where A ( t ) and B(t) are independent random processes that are autocorrelated in time through a window which is mainly associated with the sample volume shape (refer to [ 141 for a detailed description of the model). The window was an asymmetric triangle. White Gaussian noise of variance 0; was added to each of the quadrature signals. With at the variance of the simulated signal, the signal-to-noise ratio (additive SNR) is defined as SNRds = 10 log a:/a:. Several SNR’s ranging from -20 to 30 dB were generated, and power spectra (S(w)) were computed with four algorithms: FFT, Burg, MCM, and MLM. For the first three of these spectra, the Doppler frequency was computed both by a normalized first moment method, i.e., wD = SwS(w) d w / S S ( w ) dw, and by selecting the frequency of maximum spectral amplitude. The MLM spectra are not true power spectra [15], and hence it was considered inappropriate to apply the first moment to these. However, the frequency of maximum amplitude was studied for these spectra.
C. Signals Derived from Flow Experiments
The steady flow system shown in Fig. 1 was used to measure stationary Doppler signals. A constant pressure head in reservoir 1 provides flow through a 2.54-cm (l-in) diameter tube. The downstream end of the tube empties into another constant head reservoir. Fluid is returned to the upstream tank through a centrifugal pump. Flow rate is controlled by the height of the upstream tank. Test sections were inserted near the downstream end of the system. Two different sections were used: a 1-in-diameter extension of the system tube, and a straight 3.19-mm-diameter tube. The I-in test section provided a relatively B. Simulated Signals blunt velocity profile inside the sample volumes of the The effects of noise on the behavior of the spectral es- Doppler probes (velocities did not differ by more than timators were studied mainly with simulated signals. This 0.6% for the largest sample volume), whereas the has the advantage that signals with precisely known sig- 3.19-mm tube was used to study shear layers, i.e., broadband signals. The fluid was a mixture of glycerine and nal-to-noise ratios (SNR’s) can easily be generated. The model chosen for this purpose has been proposed water adjusted such that both Reynolds numbers (Re) and by Jones and Giddens [14] and generates stationary sig- velocities lay in the physiological range (viscosity was 9 nals which model the quadrature outputs obtained exper- cP). The sound speed in this fluid was 175 000 cm/s. imentally under steady laminar Row with a blunt velocity Corn starch particles with diameters ranging from 5 to 20 profile. Unlike other models proposed in the literature [ 7 ] , pm served as acoustic scatterers. To obtain parabolic flow, it incorporates not only the characteristic amplitude mod- entrance lengths were longer than L,. = (0.08Re 0.7)D ulation of the Doppler signal, but also the random phase [ 161, where D is the diameter of the tube. Flow rates were modulation which arises from the scattering process. This measured from the time required to collect a known
+
DAVID cr
U/
59 I
MODERN SPECTRAL A N A L Y S I S T k C H N I Q U k S
Drain
1I
Fig. I . Flow system employed to obtain Dopplcr data
amount of fluid. Error in this method is estimated as 1 % . The velocity at the center of the probe sample volume was measured directly with laser Doppler anemometry. For the large model, the heavy seeding required for the ultrasound measurements precluded simultaneous ultrasound and laser measurements because the laser light was too strongly attenuated by the seeding. Thus. laser measurements were taken with the probe in position but prior to ultrasonic seeding. For the small model. the shorter distance of beam travel resulted in less attenuation and measurements could be taken through the seeded fluid. Two Doppler units were used for the experimental studies: a 10-MHz system designed by C. Hartley at the Instrumentation Development Laboratories, Baylor College of Medicine, Houston, TX, and a 20-MHz Doppler catheter system manufactured by Millar Instruments, Inc. The different characteristics of the probes led to different SNR's, as will be discussed in the next section. The probes were positioned inside the tube and held by a micrometer device. The angle between probe axis and the direction of flow was zero degrees, and the sample volume was located upstream from the tip of the probe to avoid flow disturbances from the probe (refer to Fig. I ) . SNR for a given signal was computed from the ratio of the power between 0.8 and 1.2 of the peak of the Fourier spectrum divided by the power outside this band. This "in-band" measure was considered more appropriate than the ratio of total power to the power under zero flow conditions because the broad-band noise in the Doppler signals involves more than just instrument noise. For comparison with the additive noise measure defined earlier, the simulated signals were characterized by both measures. The Doppler quadrature signals Q, and Q2 were digitized on-line with a sampling rate of 30 kHz (unless stated otherwise) and stored in a Masscomp 5500 computer. Before sampling, the signals were low-pass filtered with a 13-kHz cutoff to avoid aliasing. The complex signal Q , + iQz was then analyzed off-line as described above for the simulated signals. 111. RESULTS Simulated and real signals were processed to yield velocity data as a function of time. In the data reduction,
8.53-ms records were used, which yielded 256 points per record at the sampling rate of 30 kHz. All FFT's were zero-padded to 1024 points to interpolate the spectral shape. The effects of record length and fluid velocity on the behavior of the velocity estimators are similar to those described by Jones and Giddens [ 101 in their comparison of AR and FFT methods. Since the MLM and MCM were found to behave like the Burg approach with respect to these two parameters, these effects will not be specifically addressed here. Since the signals are stationary, an ideal velocity meter would show straight horizontal lines for plots of velocity as a function of time. The performance of the algorithms was therefore assessed on the basis of two criteria: the bias and the coefficient of variation of the estimates. A . Effects @Noise
The signal-to-noise ratio of Doppler signals differs from one instrument to another for various reasons. These include differences in electronics design, signal attenuation in the medium, which is strongly frequency dependent, and probe beam patterns. The 20-MHz catheter system studied here, for example, had a much lower signal-tonoise ratio than did the 10-MHz system, as can be seen in Fig. 2. Both of the 256 point spectra in this figure are from data taken under conditions of uniform, laminar steady flow with 10- and 20-MHz probes directed upstream. Power was normalized by the total power in the signal. Broad-band noise is evident in Fig. 2(b). As the signal becomes more attenuated with range, this noise becomes more dominant. Example velocity-time traces are shown in Fig. 3 for moderate SNR. These data were from the flow system and were taken with the 20-MHz catheter system at a range of 6 mm. The flow velocity was 45.3 cm/s. The dashed line represents the time averaged velocity from the laser Doppler anemometry measurements. The FFT first moment shows strong bias and both FFT techniques show wide variance for this SNR. Similar results were obtained by Jones and Giddens [ 101 for signals with colored noise. The biases for the Burg, MLM, and FFT peak trackers were 3.1 % , 4.4 % , and 3.1 % , respectively; all were small and comparable. The MLM frequencies were consistently lower than the Burg frequencies by about 1.5% to 2 % throughout the flow data. Fig. 4 shows how the coefficient of variation of each method changes with SNR and model order. Again, the profile was blunt and flow was laminar with a velocity of 45.3 cm/s. The SNR for the 20-MHz catheter was changed via the range control. The figures are thus for decreasing SNR from top to bottom, with the highest SNR for the IO-MHz device. The coefficient of variation for the FFT is given by a horizontal line in each figure. The minimum for each curve represents the lowest coefficient of variation and hence the best model order for each algorithm. For the peak tracker, there was always a model order for the modern techniques which resulted in a lower
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 6, JUNE 1991
592
E D
0
5000.
10000.
0.
5000.
10000. (Hz)
15000
0.
Frequency
Frequency ( H r )
5.
(a) (b) Fig. 2 . Typical spectra obtained under identical experimental conditions and showing different SNR's: ( a ) IO-MHz probe: ( b ) 20-MHz catheter.
10.
15.
20.
0.
5.
10.
15.
20.
15.
20.
15.
20.
(b)
(a)
I " " I " " 1 " l
1 " " l " " I " I
t 0.
1 5.
10.
15.
20.
0.
5.
10.
(d)
(c)
0
0.
5.
10. 15. Model O r d e r
(e) 0
5 T i m e (sec)
1.
0.
.5 1 T i m e (sec)
(c) (d) Fig. 3 . Frequency-time curves for experimental signals in moderate noise: (a) Burg AR peak tracker. (b) MLM peak tracker. (c) FFT peak tracker. (d) FFT first moment. Dashed line-Doppler shift computed from LDA measurement.
coefficient of variation than the FFT peak tracker. For the high signal-to-noise case any first moment was as good as the best modem peak tracker [Fig. 4(a), (b)]. As noise increased, the precision of the first moment estimate degraded. The accuracy degraded as well. The minimum coefficient of variance for the modern peak trackers showed little dependence on SNR, but the model order associated with this minimum increased with decreasing SNR. An order of 5 performed well for all SNR's studied. In this respect, both AR methods behaved nearly identically. Where the MLM method differed from the AR methods, e.g., for the higher model orders, its variance was lower. Also, its coefficients varied more smoothly than those of the AR methods. First moments for all techniques were nearly identical except for very low model orders and noisy signals [Fig. 4(d), (f)]. A more extensive study was undertaken with the simulated signals. The phenomena described above for actual
20.
0.
5.
10.
Model O r d e r
(f)
Fig. 4. Coefficients of variation as a function of model order for three experimentally observed SNR's. (a), (c). and (e): peak trackers. (b). (d). and ( f ) : first moments (a). (b) IO-MHz probe, SNR = 21 dB. (c). (d) 20-MHz catheter, range 3 mm, SNR = 1 I dB. (e), ( f ) 20-MHz catheter. range I O mm,SNR = 9 dB.
Doppler signals agreed with those observed with the simulations. This confirmed that the effects were caused directly by the change in SNR and not by other signal characteristics which might change between two different Doppler ultrasound instruments or one instrument at two different ranges. The mean and the coefficient of variation of the estimates are shown as a function of the additive SNR in Fig. 5 . The scale above the plots shows the inband signal-to-noise measure. An order of 5 was used for the modern methods. Results for the MCM are identical to those of the Burg method and are therefore not shown. Also, first moment frequencies derived from the AR or FFT spectra were not significantly different. The timeaveraged frequencies for all of the peak trackers agreed with one another throughout the SNR range studied. For SNR's above - 10 dB, they were all near the input Doppler shift of the simulation (the horizontal line in the figure). Below this SNR, they became progressively diminished. The 1.5% to 2 % difference between the AR and MLM methods discussed in connection with the flow data
DAVID
"
-5.
I
593
MODERN SPECTRAL ANALYSIS I t . ( HNIQUES
CI U /
0. 5 . 10. 1 5 .
I
I
I
I
-5.
18.
1
\
0. 5
/
I
10. 15.
l
18
peak trackers over the first moments in the presence of colored noise.
I
l
0
t
-20.
-10.
0. 10 20 A d d i t i v e SNR ( d b )
(a)
30.
o
n
i
L
-20
-10
0 10 20 A d d i t i v e SNR ( d b )
30
(b)
Fig. 5 . Estimated ( a ) mean frequcncq and ( h ) coetficienth t i l variancc as a function of S N R for simulated signals. Model orders are :dl 5 . The lower .r-axis scale is the additive S N K , and the upper i-axis scale I S the in-hand SNR. Solid line--exact frequency : circle-til-st iiioiiienth: square-Burg peak tracker; trianSlc-MLM peak tracker: diaiiiond-FFT peak tracker.
was not seen in the simulation results. The time-averaged frequency for the first moment was significantly lower than the simulation frequency until the SNR reached 15 dB. The coefficients of variation for all the peak trackers increased dramatically with decreased SNR below 0 dB. Peak frequency variances for the AR and MLM methods were lower than those for the FFT. The first moment variances were higher than the AR and MLM peak frequency variances until SNR reached 18 dB. Of note is the increase in the variance for the Burg peak tracker at a SNR above 18 dB. For these high signal-to-noise ratios. 5 is not the optimal model order for the Burg spectrum, but rather a lower model order is preferable. This has been discussed by Jones and Giddens [ 141. The MLM is much less sensitive to model order in this respect. Results indicate that significant improvements over the FFT can be achieved with modern estimation techniques coupled with a peak tracking algorithm in moderate SNR environments. Although these improvements are greatest when the optimal filter length is used, they are still significant if a constant model order ( 5 , for example) is chosen. Inspection of spectra to determine the noise level could help in choosing the optimal processing order, but this is not practical in routine (especially medical) applications. It is therefore of interest to find basic relations between the SNR and the optimum filter length. The model selected for the noise i n the simulations may not correspond to experimental conditions. The bias in the first moment for this particular model is readily explained because white noise contains equally distributed spectral content on either side of zero frequency. Thus, the bias will always be toward zero [ IO]. For other types of noise, bias will still be present and will be even greater if the spectral content of the noise is concentrated within a small frequency band. Thus, the results obtained up to now should be valid in every low SNR environment. Several measurements taken with the IO-MHz probe positioned outside of the fluid at a 60" angle from its direction of motion revealed colored low frequency noise. The processing of these signals confirmed the superiority of the
B. Effects oj'Sigizul Bandwidth Catheter measurements were taken under shear layer conditions to determine the effect of signal bandwidth on the behavior of the velocity estimators. The probe was mounted (directed upstream) in the 3.19-mm-diameter tube and measurements were taken at 5 positions across the diameter. The data in this case were digitized at 40 kHz. Since a parabolic velocity profile was present, the sample volume always contained a range of velocities. The sample volume dimensions were measured with a string target and found to be roughly 1 mm in depth by 1 mm in diameter based on the 6-dB down location. The estimated mean velocity was approximately the same for all the first moment estimators on one hand, and for all the peak trackers on the other hand, if model orders of 1 to 3 are excluded. Mean estimated velocities and coefficients of variance for model order 5 are shown as a function of location (hence as a function of the velocity gradient) in Fig. 6. For comparison, laser Doppler measurements are included in this graph. The large probe sample volume tends to depress the estimated velocities for both ultrasound and laser measurements, particularly near the centerline. The largest peak tracker error is for the MLM methods and is 1 1 . 7 % . The first moment approach, is accurate only near the wall. This is probably independent of the presence of shear and is simply related to the lower velocity in the sample volume: the error in the first moment estimation increases with the velocity in the presence of low frequency noise. Fig. 6(b) shows that the standard deviations of the modern peak trackers with model order 5 are generally smaller than those of the FFT tracker. Model order 5 is a compromise between a small standard deviation and an accurate estimate of the average (smaller orders underestimate the velocity, whereas higher orders are subject to higher fluctuations). The data at + 1 mm in Fig. 6(b) show that for peak tracking, the variation in the order 5 Burg AR method can be comparable to that for the FFT in strong shear layers. The asymmetry visible in Fig. 6(b) was consistent over 4 independent profile measurements and may indicate an asymmetry in the probe sample volume. The results in weak or moderate shear layers (the 3 positions away from the wall. where the spectrum is narrower) are similar to those described earlier for moderate SNR's. For the strong shear layer region, velocity-time curves are shown in Fig. 7. In this sample of curves, the MLM provides the lowest variance of all the peak trackers. The FFT peak tracker frequency varies widely. The AR techniques show less variation, but occasional spikes are present which are comparable to the FFT peak tracker fluctuations. The FFT first moment variance is the lowest, but this advantage must be weighed against the noise bias of this technique. The jumps from one frequency to another for the AR peak trackers [Fig. 7(a)] are similar to the behavior of
594
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38. NO. 6, JUNE 1991 -
-
0
" m
0 N
.E 9 3
x
0
*
n
-
6
O U ;
0
U
-
~~~
0
, a U W
N 0
0
0
m
0
o
-1.
0. 1. R o d i a l P o s i t i o n (m)
0.
-1
Rodiol P o s i t i o n
iaJ
0
0.
5000.
10000.
15000:
0.
5000.
10000.
15000.
0.
5000.
10000.
15000,
0.
5000.
10000,
15000.
5000.
10000,
15000.
0.
5000.
10000.
15000
1.
(mn)
(b) -n
Fig. 6. Velocity estimation in a 3.19." diameter tube. Dot-laser Doppler measurements; circle-first moments: square-Burg peak tracker: triangle-MLM peak tracker: diamond-FFT peak tracker. (a) Time averaged velocities. (Solid line is the parabolic profile.) (b) Standard deviations. (Solid lines are to improve readability.)
0 N
;. 9 3
f 0
-
L 9
-
0
'
6 0
0.
5
1. 1.5 T i m e (sec)
2.
0.
5
1. 1 5 Time (sec)
2.
(b)
(a)
0.
F r e q u m c y (Hz)
(e)
Frequency
(Hz)
(f)
Fig. 8. Singlc (left-hand side) and averaged (right-hand side) spectral estimates. (a), (b) FFT (rectangular window): (c). (d) AR order 16 (fpe); (e), i f ) MLM order 25.
C. Spectrum Estimation 0.
.5
1.
1.5
Time (sec)
2.
0.
.5
1. 1.5 Time (sec)
2.
(C) (d) Fig. 7 . Frequency-time curves for a shear layer (a) Burg AR peak tracker. ( b ) MLM peak tracker. ( c ) FFT peak tracker, id) FFT first moment.
phase locked loop systems [ 171. For such broad-band signals, a tracking algorithm can be used which takes into account the past history of the velocity signal and disallows sudden changes in the Doppler frequency. However, such algorithms are sensitive in that one looses the track easily when the velocity changes rapidly, a situation encountered physiologically. The MLM provides the most stable peak tracker. Although it is associated with a small underestimation of the fluid velocity, it is nearly as stable as the first moment estimators. Moreover, AR peak trackers have been found to be very sensitive to the choice of model order, even in the range 3-7, whereas the MLM tracker was more robust. Thus, the MLM peak tracker is preferable if the nature of signal (broad or narrow band) is not known.
It is desirable to display the spectrum of the signal because clinicians take a qualitative impression of spectral broadening into account during diagnosis [ 181. Under any experimental conditions, but especially in a shear layer, inspection of single spectra reveals large fluctuations in the spectral shape with time, even though the flow is steady. Spectral averaging greatly reduces these fluctuations and yields spectra which are consistent with the velocity profile independently of the spectral estimation method. However, this technique is difficult to implement reliably in pulsatile flow. The question thus arises as to which algorithm yields the most stable spectral estimates. Fig. 8 shows averaged and single spectra for data recorded in the strong shear layer. The spectra on the right [Fig. 8(b), (d), and (f)] are ensembles over 64 records. The figures on the left [Fig. 8(a), (c), and (e)] are 12 individual spectra plotted over one another. The three ensembled spectra have similar shapes. The unensembled realizations of the FFT spectrum in Fig. 8(a) show the widest variation. The AR spectra in Fig. 8(c) show less variability, and the MLM spectra in Fig. 8(e) show the least variability.
D A V I D er
< I / , :
595
MODERN SPECTRAL ANALYSIS TECHNIQUES
The MCM behave similarly to the Burg algorithm for both single and averaged spectra, independently of the model order, and hence the MCM spectra are not displayed. Furthermore, low orders for the modem methods result in unrealistic spectral shapes. Windowing the FFT does not significantly change either single or averaged spectra. For the AR methods, all orders above the FPE yield results similar to the one plotted. The FPE indicated orders between 12 and 18. For the MLM, an order of 25 was ultimately chosen because it yielded spectra consistent with the velocity distribution inside the sample volume for every degree of shear encountered. Higher orders were statistically less stable and computationally less efficient. IV. DISCUSSION The AR and MLM techniques have been shown to be superior to Fourier methods under nearly all experimental conditions when a peak tracking algorithm is used to find the Doppler frequency, provided the model order is chosen appropriately. An order of 5 performed well under all flow, noise, and signal bandwidth conditions investigated. Applying a Hanning window to the signal did not improve the behavior of the FFT. The performance of all the peak trackers was dependent on the bandwidth of the signal, i.e., the range of velocities present in the sample volume. The MLM yielded the most robust peak tracker in that it was the least sensitive to changes in model order and signal bandwidth. The Burg algorithm appears to be more robust than the MCM for peak tracking with the Doppler signals. For model orders above 15, numerical ill conditioning sometimes occurred for the MCM, and the MCM would also result in peaks which did not represent Doppler frequencies. These peaks were sometimes stronger than the actual Doppler shift peaks and led to spikes in the velocity-time curves from the peak detector. When the first moment technique was used to compute the velocity, the autoregressive and Fourier methods were virtually identical. The computation was nearly independent of the model order for the AR algorithms, and again, nonrectangular window shapes did not improve the FFT results. The first moment approach was barely sensitive to the signal bandwidth, but strongly affected by the presence of noise. The relative performance of the peak tracking and first moment methods depends on the level of noise in the signal and on the range of velocities present in the sample volume. The peak trackers were much more accurate (both in terms of bias and variance) than the first moments when noise was present in the signal. This was observed for both colored and white noise. However, they became unstable when the spectrum of the signal was broad, unlike the first moments. For spectral estimation, the MLM algorithm appeared to be more stable than the other approaches. Its use should be particularly rewarding in the analysis of pulsatile flow
signals, where the estimation of the velocity distribution inside the sample volume relies on single spectral estimates because reliable spectral averaging is difficult. It must be remembered, however, that the MLM spectrum is not a true power spectrum. The area under the spectrum between two frequencies is not the power in that frequency band. However, the shape of the spectrum still provides insight into the velocity field within the sample volume. Care must be taken if these spectra are to be used to extract specific power related indices. Comparison of the results presented here to those of Vaitkus and Cobbold [9] must be made cautiously because different model orders and different algorithms were used in the computations. Also, the spectra here apply to pulsed Doppler ultrasound whereas those of Vaitkus and Cobbold model continuous wave spectra. However, the lower spectral variance for the MLM method over both Fourier and AR methods is notable. The difference in optimal AR model order for the two studies, 16 for this study, 8 for Vaitkus and Cobbold reflects the broader nature of the spectra studied here. Questions still exist in regard to the applicability of the modern spectral analysis techniques to more complicated spectra, such as might be obtained from a stenosed artery. The following strategy for processing pulsed Doppler ultrasound signals can be proposed from the above observations: display simultaneously the MLM order 5 peak tracker, a first moment estimate, and a higher order MLM spectrum (an order of 25 seemed to be satisfactory for all the flow conditions studied). This procedure allows an estimate of the bandwidth of the signal, and thus a choice of the best suited velocity estimator (the peak tracker should always be preferred except for very broad spectra). This strategy is being proposed with currently available Doppler instruments equipped with FFT processors [ 191. It is clear that such a strategy can be optimized by using modem spectral analysis algorithms instead of FFT’s, since the AR yields first moments identical to those of the FFT, while the MLM peak tracker is always superior to the FFT tracker when the order is kept below IO. Since the MLM spectrum is computed from the AR parameters, it is not unrealistic to compute the first moment of the Burg spectrum while displaying the MLM spectrum. V. CONCLUSION Two important parameters influencing the accuracy of velocity estimation from pulsed Doppler ultrasound signals as obtained with four spectral analysis algorithms have been studied. It was shown that the optimal algorithm is dependent on the flow (bandwidth) and the signalto-noise conditions. Robustness considerations lead to the choice of a combined display of MLM order 5 peak tracker and high order (25) MLM spectrum. An AR first moment estimator can also be displayed if broad-band signals are expected. The results may be clinically significant since accurate and robust velocity estimators are needed in order to de-
596
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 38. NO. 6. J U N E 1991
rive reliable diagnostic ~ l visual~ inspection ~ , of spectra as well as quantitative studies on spectral shape (e.g., spectral broadening index) require stable and accurate spectral estimators.
REFERENCES W . K . George and J . J . Luniley. "The laser Doppler velocimeter and its application to the measurement of turbulence," J . Fl~tidM c d i . . vol. 60, pp. 321-362. 1973. S. M. Kay and S. L. Marple, "Spectrum analysis--a modern perspective." Proc. /EEL-. vol. 69. pp. 1380-1419. 1981. L. J . D'Luna and V . L. Newhouse. "Vortex characterization and
1171 D. P. Giddens and A . M. A . Khalifa, "Turbulence nieasureinents u'ith pulsed Doppler ultrasound employing a frequency tracking method.'' Ultrctsoir~ltlM o d . Biol.. vol. 8. pp. 427-437. 1982. 1181 W . J . Zwiebel. "Spectrum analysis in carotid sonography." Ultrci,sourid M r d . Biol.. vol. 13. pp. 625-636. 1987. 1191 W . F . Voyles. S. A . Altobclli, D. C. Fisher, and E. R. Greene. " A comparison o l digital and analog inethods of Doppler spectral analysis for quantifying flow." Ulrro.sorirrtl M e d . Biol.. vol. I I pp. 727733. 10x5. I
Jean-Yves David received the engineering degree from the University of Technology of Comidentification by ultrasound Doppler," Ulrrct.soiiic. /i~~ei
589
Modern Spectral Analysis Techniques for Blood Flow Velocity and Spectral Measurements with Pulsed Doppler Ultrasound Jean-Yves David, Steven A. Jones, Member, IEEE, and Don P. Giddens
Abstract-Four spectral analysis techniques were applied to pulsed Doppler ultrasonic quadrature signals to compare the relative merits of each technique for estimation of flow velocity and Doppler spectra. The four techniques were 1) the fast Fourier transform method, 2) the maximum likelihood method, 3) the Burg autoregressive algorithm, and 4) the modified covariance approach to autoregressive modeling. Both simulated signals and signals obtained from an in vitro flow system were studied. Optimal parameter values (e.g., model orders) were determined for each method, and the effects of signal-to-noise ratio and signal bandwidth were investigated. The modern spectral analysis techniques were shown to be superior to Fourier techniques in most circumstances, provided the model order was chosen appropriately. Robustness considerations tended to recommend the maximum likelihood method for both velocity and spectral estimation. Despite the restrictions of steady laminar flow,, the results provide important basic information concerning the applicability of modern spectral analysis techniques to Doppler ultrasonic evaluation of arterial disease.
I . INTRODUCTION
P
ULSED Doppler ultrasound instruments have been used in clinical applications which range from the study of fetal blood flow to the diagnosis of atherosclerosis. In one of the most common signal analysis methods used with these instruments, the received pulses are downmixed to quadrature signals and these are transformed into power spectra. Although the principle of Doppler measurement is simple, the spectral content of quadrature signals is complicated by a number of devicerelated factors which introduce uncertainty in the velocity estimations. Low spatial resolution in regions of large velocity gradients, low temporal resolution for rapidly changing flows, instrument noise, and transit time effects all contribute to a broadening of the Doppler power spectra [ I ] . However, such broadening may also arise from limitations in the resolution and accuracy of the spectrum estimator employed. It is thus important to determine the appropriateness of the available spectrum estimator models to Doppler ultrasound signals. Manuscript received December 7. 1989: rebised August 13, 1990. J:Y. David is with the School o f Mechanical Engineering, Georgia Institute o f Technology. Atlanta. C A 30332. S . A. Jones and D. P. Giddens are with thc School o l Aerospace Engineering, Georgia Institute of Technology, Atlanta. C A 30332. IEEE Log Number 9144684.
A number of spectral estimation techniques have recently been developed and have been compared to the more standard fast Fourier transform (FFT) method for Doppler ultrasonic signal processing. The most common of these methods are autoregressive (AR), moving average, autoregressive moving average, and the maximum likelihood method (MLM). In addition, AR spectra can be computed via several algorithms. These include the Burg method, the Yule-Walker method, and the modified covariance method (MCM). These techniques are described in more detail by Kay and Marple [2]. D’luna et al. [3], [4] used the Burg method for AR modeling to investigate in virro poststenotic vortices. The final prediction error criterion (FPE) of Akaike [ 5 ] was used to determine the appropriate model order. Use of the AR was motivated by its greater spectral resolution, but no comparison of the autoregressive approach with other methods was attempted. Kaluzynski [6] compared FFT and AR methods and showed that for short data records, the AR yields more statistical stability than the FFT. He did not address the accuracy of the velocity estimation, but considered only the stability of the spectral shape. Again, the FPE was used for the model order selection. Talhami and Kitney [7] investigated an autoregressive moving average model but did not address the problem of model order selection. Moreover, they used simulated signals that did not show the random phase fluctuations characteristic of real Doppler signals. An extensive study of several modem spectral analysis methods, including AR, MLM, and autoregressive moving average, was done by Vaitkus and Cobbold [8], [9]. Simulated signals were used which modeled the shape of a spectrum taken in vivo with continuous wave ultrasound, but the simulation was not based on the physics of the Doppler process. The study was directed toward spectral estimation and did not discuss velocity estimation. Limited in vivo data were presented, and although the actual in vivo flow field was unknown, many of the results of the simulation studies were supported by these data. For velocity estimation, Jones and Giddens [ 101 provided a comparison between FFT methods and Burg AR modeling. This work can be viewed as an extension of these investigations. Few of the investigations to date have quantitatively examined the accuracy of the modern spectral analysis
0018-929419 110600-0589$01.OO O 1991 IEEE
590
I E E E T R A N S A C T I O N S O N BlOMEDlCAL E N G I N E E R I N G , VOL. 38. N O . 6, J U N E 1991
techniques. Those investigations which were quantitative have been undertaken with simulations rather than actual flow data. These simulations were designed to model spectral shape or amplitude modulation phenomena rather than the underlying physics of Doppler signals. Although much work has been done to examine spectral estimation, velocity estimation methods with modern techniques have not been extensively studied. The effect of signal-to-noise ratio on the accuracy of velocity measurements has also been neglected. This paper compares the performance of the FFT, MLM, MCM, and Burg AR under a variety of signall noise and signal bandwidth conditions. Simulated signals which incorporate the essential physics of Doppler measurement and signals derived from laboratory flow studies will be used for this purpose. Both spectral and velocity estimates will be examined. The spectral analysis methods to be studied were selected on the following basis. The Burg AR technique is known to perform well with signals of short duration; MCM exhibits globally better behavior than Burg AR for noisy sinusoids; MLM has a smaller variance than the autoregressive methods under some conditions [ 1 11, [ 121. The FFT was used for comparison since it is the technique commonly employed in commercially available devices. 11. METHODS A . Spectral Analysis Methods The spectral analysis methods employed have been described extensively in the literature (21, [121. The periodogram was used with either a rectangular or a Hanning window. Autoregressive spectra were obtained via both the Burg algorithm and the MCM approach. MLM spectra were evaluated by the method of Musicus [ 131, which uses the AR parameters from the Burg algorithm to estimate the spectrum. Each algorithm is implemented for complex signals in order to retain the directionality information of the Doppler quadrature signals. The convergence criteria for the MCM were switched off during the following comparisons so that the program could iterate until a specified order had been reached or until it terminated automatically upon detection of numerical difficulty. A detailed description of each estimator may be found in 121 and 112).
modulation is manifested by random fluctuations in the detected Doppler shifted frequency that are not related to actual flow fluctuations. This phenomenon has been termed “Doppler ambiguity” [ l ] , and it affects the velocity and spectral estimations. Briefly, the simulation model generates two quadrature signals expressed as
Q l ( t ) = A(t) cos
(U(,[)
+ B ( f )sin (w,,t)
(1)
where A ( t ) and B(t) are independent random processes that are autocorrelated in time through a window which is mainly associated with the sample volume shape (refer to [ 141 for a detailed description of the model). The window was an asymmetric triangle. White Gaussian noise of variance 0; was added to each of the quadrature signals. With at the variance of the simulated signal, the signal-to-noise ratio (additive SNR) is defined as SNRds = 10 log a:/a:. Several SNR’s ranging from -20 to 30 dB were generated, and power spectra (S(w)) were computed with four algorithms: FFT, Burg, MCM, and MLM. For the first three of these spectra, the Doppler frequency was computed both by a normalized first moment method, i.e., wD = SwS(w) d w / S S ( w ) dw, and by selecting the frequency of maximum spectral amplitude. The MLM spectra are not true power spectra [15], and hence it was considered inappropriate to apply the first moment to these. However, the frequency of maximum amplitude was studied for these spectra.
C. Signals Derived from Flow Experiments
The steady flow system shown in Fig. 1 was used to measure stationary Doppler signals. A constant pressure head in reservoir 1 provides flow through a 2.54-cm (l-in) diameter tube. The downstream end of the tube empties into another constant head reservoir. Fluid is returned to the upstream tank through a centrifugal pump. Flow rate is controlled by the height of the upstream tank. Test sections were inserted near the downstream end of the system. Two different sections were used: a 1-in-diameter extension of the system tube, and a straight 3.19-mm-diameter tube. The I-in test section provided a relatively B. Simulated Signals blunt velocity profile inside the sample volumes of the The effects of noise on the behavior of the spectral es- Doppler probes (velocities did not differ by more than timators were studied mainly with simulated signals. This 0.6% for the largest sample volume), whereas the has the advantage that signals with precisely known sig- 3.19-mm tube was used to study shear layers, i.e., broadband signals. The fluid was a mixture of glycerine and nal-to-noise ratios (SNR’s) can easily be generated. The model chosen for this purpose has been proposed water adjusted such that both Reynolds numbers (Re) and by Jones and Giddens [14] and generates stationary sig- velocities lay in the physiological range (viscosity was 9 nals which model the quadrature outputs obtained exper- cP). The sound speed in this fluid was 175 000 cm/s. imentally under steady laminar Row with a blunt velocity Corn starch particles with diameters ranging from 5 to 20 profile. Unlike other models proposed in the literature [ 7 ] , pm served as acoustic scatterers. To obtain parabolic flow, it incorporates not only the characteristic amplitude mod- entrance lengths were longer than L,. = (0.08Re 0.7)D ulation of the Doppler signal, but also the random phase [ 161, where D is the diameter of the tube. Flow rates were modulation which arises from the scattering process. This measured from the time required to collect a known
+
DAVID cr
U/
59 I
MODERN SPECTRAL A N A L Y S I S T k C H N I Q U k S
Drain
1I
Fig. I . Flow system employed to obtain Dopplcr data
amount of fluid. Error in this method is estimated as 1 % . The velocity at the center of the probe sample volume was measured directly with laser Doppler anemometry. For the large model, the heavy seeding required for the ultrasound measurements precluded simultaneous ultrasound and laser measurements because the laser light was too strongly attenuated by the seeding. Thus. laser measurements were taken with the probe in position but prior to ultrasonic seeding. For the small model. the shorter distance of beam travel resulted in less attenuation and measurements could be taken through the seeded fluid. Two Doppler units were used for the experimental studies: a 10-MHz system designed by C. Hartley at the Instrumentation Development Laboratories, Baylor College of Medicine, Houston, TX, and a 20-MHz Doppler catheter system manufactured by Millar Instruments, Inc. The different characteristics of the probes led to different SNR's, as will be discussed in the next section. The probes were positioned inside the tube and held by a micrometer device. The angle between probe axis and the direction of flow was zero degrees, and the sample volume was located upstream from the tip of the probe to avoid flow disturbances from the probe (refer to Fig. I ) . SNR for a given signal was computed from the ratio of the power between 0.8 and 1.2 of the peak of the Fourier spectrum divided by the power outside this band. This "in-band" measure was considered more appropriate than the ratio of total power to the power under zero flow conditions because the broad-band noise in the Doppler signals involves more than just instrument noise. For comparison with the additive noise measure defined earlier, the simulated signals were characterized by both measures. The Doppler quadrature signals Q, and Q2 were digitized on-line with a sampling rate of 30 kHz (unless stated otherwise) and stored in a Masscomp 5500 computer. Before sampling, the signals were low-pass filtered with a 13-kHz cutoff to avoid aliasing. The complex signal Q , + iQz was then analyzed off-line as described above for the simulated signals. 111. RESULTS Simulated and real signals were processed to yield velocity data as a function of time. In the data reduction,
8.53-ms records were used, which yielded 256 points per record at the sampling rate of 30 kHz. All FFT's were zero-padded to 1024 points to interpolate the spectral shape. The effects of record length and fluid velocity on the behavior of the velocity estimators are similar to those described by Jones and Giddens [ 101 in their comparison of AR and FFT methods. Since the MLM and MCM were found to behave like the Burg approach with respect to these two parameters, these effects will not be specifically addressed here. Since the signals are stationary, an ideal velocity meter would show straight horizontal lines for plots of velocity as a function of time. The performance of the algorithms was therefore assessed on the basis of two criteria: the bias and the coefficient of variation of the estimates. A . Effects @Noise
The signal-to-noise ratio of Doppler signals differs from one instrument to another for various reasons. These include differences in electronics design, signal attenuation in the medium, which is strongly frequency dependent, and probe beam patterns. The 20-MHz catheter system studied here, for example, had a much lower signal-tonoise ratio than did the 10-MHz system, as can be seen in Fig. 2. Both of the 256 point spectra in this figure are from data taken under conditions of uniform, laminar steady flow with 10- and 20-MHz probes directed upstream. Power was normalized by the total power in the signal. Broad-band noise is evident in Fig. 2(b). As the signal becomes more attenuated with range, this noise becomes more dominant. Example velocity-time traces are shown in Fig. 3 for moderate SNR. These data were from the flow system and were taken with the 20-MHz catheter system at a range of 6 mm. The flow velocity was 45.3 cm/s. The dashed line represents the time averaged velocity from the laser Doppler anemometry measurements. The FFT first moment shows strong bias and both FFT techniques show wide variance for this SNR. Similar results were obtained by Jones and Giddens [ 101 for signals with colored noise. The biases for the Burg, MLM, and FFT peak trackers were 3.1 % , 4.4 % , and 3.1 % , respectively; all were small and comparable. The MLM frequencies were consistently lower than the Burg frequencies by about 1.5% to 2 % throughout the flow data. Fig. 4 shows how the coefficient of variation of each method changes with SNR and model order. Again, the profile was blunt and flow was laminar with a velocity of 45.3 cm/s. The SNR for the 20-MHz catheter was changed via the range control. The figures are thus for decreasing SNR from top to bottom, with the highest SNR for the IO-MHz device. The coefficient of variation for the FFT is given by a horizontal line in each figure. The minimum for each curve represents the lowest coefficient of variation and hence the best model order for each algorithm. For the peak tracker, there was always a model order for the modern techniques which resulted in a lower
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 6, JUNE 1991
592
E D
0
5000.
10000.
0.
5000.
10000. (Hz)
15000
0.
Frequency
Frequency ( H r )
5.
(a) (b) Fig. 2 . Typical spectra obtained under identical experimental conditions and showing different SNR's: ( a ) IO-MHz probe: ( b ) 20-MHz catheter.
10.
15.
20.
0.
5.
10.
15.
20.
15.
20.
15.
20.
(b)
(a)
I " " I " " 1 " l
1 " " l " " I " I
t 0.
1 5.
10.
15.
20.
0.
5.
10.
(d)
(c)
0
0.
5.
10. 15. Model O r d e r
(e) 0
5 T i m e (sec)
1.
0.
.5 1 T i m e (sec)
(c) (d) Fig. 3 . Frequency-time curves for experimental signals in moderate noise: (a) Burg AR peak tracker. (b) MLM peak tracker. (c) FFT peak tracker. (d) FFT first moment. Dashed line-Doppler shift computed from LDA measurement.
coefficient of variation than the FFT peak tracker. For the high signal-to-noise case any first moment was as good as the best modem peak tracker [Fig. 4(a), (b)]. As noise increased, the precision of the first moment estimate degraded. The accuracy degraded as well. The minimum coefficient of variance for the modern peak trackers showed little dependence on SNR, but the model order associated with this minimum increased with decreasing SNR. An order of 5 performed well for all SNR's studied. In this respect, both AR methods behaved nearly identically. Where the MLM method differed from the AR methods, e.g., for the higher model orders, its variance was lower. Also, its coefficients varied more smoothly than those of the AR methods. First moments for all techniques were nearly identical except for very low model orders and noisy signals [Fig. 4(d), (f)]. A more extensive study was undertaken with the simulated signals. The phenomena described above for actual
20.
0.
5.
10.
Model O r d e r
(f)
Fig. 4. Coefficients of variation as a function of model order for three experimentally observed SNR's. (a), (c). and (e): peak trackers. (b). (d). and ( f ) : first moments (a). (b) IO-MHz probe, SNR = 21 dB. (c). (d) 20-MHz catheter, range 3 mm, SNR = 1 I dB. (e), ( f ) 20-MHz catheter. range I O mm,SNR = 9 dB.
Doppler signals agreed with those observed with the simulations. This confirmed that the effects were caused directly by the change in SNR and not by other signal characteristics which might change between two different Doppler ultrasound instruments or one instrument at two different ranges. The mean and the coefficient of variation of the estimates are shown as a function of the additive SNR in Fig. 5 . The scale above the plots shows the inband signal-to-noise measure. An order of 5 was used for the modern methods. Results for the MCM are identical to those of the Burg method and are therefore not shown. Also, first moment frequencies derived from the AR or FFT spectra were not significantly different. The timeaveraged frequencies for all of the peak trackers agreed with one another throughout the SNR range studied. For SNR's above - 10 dB, they were all near the input Doppler shift of the simulation (the horizontal line in the figure). Below this SNR, they became progressively diminished. The 1.5% to 2 % difference between the AR and MLM methods discussed in connection with the flow data
DAVID
"
-5.
I
593
MODERN SPECTRAL ANALYSIS I t . ( HNIQUES
CI U /
0. 5 . 10. 1 5 .
I
I
I
I
-5.
18.
1
\
0. 5
/
I
10. 15.
l
18
peak trackers over the first moments in the presence of colored noise.
I
l
0
t
-20.
-10.
0. 10 20 A d d i t i v e SNR ( d b )
(a)
30.
o
n
i
L
-20
-10
0 10 20 A d d i t i v e SNR ( d b )
30
(b)
Fig. 5 . Estimated ( a ) mean frequcncq and ( h ) coetficienth t i l variancc as a function of S N R for simulated signals. Model orders are :dl 5 . The lower .r-axis scale is the additive S N K , and the upper i-axis scale I S the in-hand SNR. Solid line--exact frequency : circle-til-st iiioiiienth: square-Burg peak tracker; trianSlc-MLM peak tracker: diaiiiond-FFT peak tracker.
was not seen in the simulation results. The time-averaged frequency for the first moment was significantly lower than the simulation frequency until the SNR reached 15 dB. The coefficients of variation for all the peak trackers increased dramatically with decreased SNR below 0 dB. Peak frequency variances for the AR and MLM methods were lower than those for the FFT. The first moment variances were higher than the AR and MLM peak frequency variances until SNR reached 18 dB. Of note is the increase in the variance for the Burg peak tracker at a SNR above 18 dB. For these high signal-to-noise ratios. 5 is not the optimal model order for the Burg spectrum, but rather a lower model order is preferable. This has been discussed by Jones and Giddens [ 141. The MLM is much less sensitive to model order in this respect. Results indicate that significant improvements over the FFT can be achieved with modern estimation techniques coupled with a peak tracking algorithm in moderate SNR environments. Although these improvements are greatest when the optimal filter length is used, they are still significant if a constant model order ( 5 , for example) is chosen. Inspection of spectra to determine the noise level could help in choosing the optimal processing order, but this is not practical in routine (especially medical) applications. It is therefore of interest to find basic relations between the SNR and the optimum filter length. The model selected for the noise i n the simulations may not correspond to experimental conditions. The bias in the first moment for this particular model is readily explained because white noise contains equally distributed spectral content on either side of zero frequency. Thus, the bias will always be toward zero [ IO]. For other types of noise, bias will still be present and will be even greater if the spectral content of the noise is concentrated within a small frequency band. Thus, the results obtained up to now should be valid in every low SNR environment. Several measurements taken with the IO-MHz probe positioned outside of the fluid at a 60" angle from its direction of motion revealed colored low frequency noise. The processing of these signals confirmed the superiority of the
B. Effects oj'Sigizul Bandwidth Catheter measurements were taken under shear layer conditions to determine the effect of signal bandwidth on the behavior of the velocity estimators. The probe was mounted (directed upstream) in the 3.19-mm-diameter tube and measurements were taken at 5 positions across the diameter. The data in this case were digitized at 40 kHz. Since a parabolic velocity profile was present, the sample volume always contained a range of velocities. The sample volume dimensions were measured with a string target and found to be roughly 1 mm in depth by 1 mm in diameter based on the 6-dB down location. The estimated mean velocity was approximately the same for all the first moment estimators on one hand, and for all the peak trackers on the other hand, if model orders of 1 to 3 are excluded. Mean estimated velocities and coefficients of variance for model order 5 are shown as a function of location (hence as a function of the velocity gradient) in Fig. 6. For comparison, laser Doppler measurements are included in this graph. The large probe sample volume tends to depress the estimated velocities for both ultrasound and laser measurements, particularly near the centerline. The largest peak tracker error is for the MLM methods and is 1 1 . 7 % . The first moment approach, is accurate only near the wall. This is probably independent of the presence of shear and is simply related to the lower velocity in the sample volume: the error in the first moment estimation increases with the velocity in the presence of low frequency noise. Fig. 6(b) shows that the standard deviations of the modern peak trackers with model order 5 are generally smaller than those of the FFT tracker. Model order 5 is a compromise between a small standard deviation and an accurate estimate of the average (smaller orders underestimate the velocity, whereas higher orders are subject to higher fluctuations). The data at + 1 mm in Fig. 6(b) show that for peak tracking, the variation in the order 5 Burg AR method can be comparable to that for the FFT in strong shear layers. The asymmetry visible in Fig. 6(b) was consistent over 4 independent profile measurements and may indicate an asymmetry in the probe sample volume. The results in weak or moderate shear layers (the 3 positions away from the wall. where the spectrum is narrower) are similar to those described earlier for moderate SNR's. For the strong shear layer region, velocity-time curves are shown in Fig. 7. In this sample of curves, the MLM provides the lowest variance of all the peak trackers. The FFT peak tracker frequency varies widely. The AR techniques show less variation, but occasional spikes are present which are comparable to the FFT peak tracker fluctuations. The FFT first moment variance is the lowest, but this advantage must be weighed against the noise bias of this technique. The jumps from one frequency to another for the AR peak trackers [Fig. 7(a)] are similar to the behavior of
594
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38. NO. 6, JUNE 1991 -
-
0
" m
0 N
.E 9 3
x
0
*
n
-
6
O U ;
0
U
-
~~~
0
, a U W
N 0
0
0
m
0
o
-1.
0. 1. R o d i a l P o s i t i o n (m)
0.
-1
Rodiol P o s i t i o n
iaJ
0
0.
5000.
10000.
15000:
0.
5000.
10000.
15000.
0.
5000.
10000.
15000,
0.
5000.
10000,
15000.
5000.
10000,
15000.
0.
5000.
10000.
15000
1.
(mn)
(b) -n
Fig. 6. Velocity estimation in a 3.19." diameter tube. Dot-laser Doppler measurements; circle-first moments: square-Burg peak tracker: triangle-MLM peak tracker: diamond-FFT peak tracker. (a) Time averaged velocities. (Solid line is the parabolic profile.) (b) Standard deviations. (Solid lines are to improve readability.)
0 N
;. 9 3
f 0
-
L 9
-
0
'
6 0
0.
5
1. 1.5 T i m e (sec)
2.
0.
5
1. 1 5 Time (sec)
2.
(b)
(a)
0.
F r e q u m c y (Hz)
(e)
Frequency
(Hz)
(f)
Fig. 8. Singlc (left-hand side) and averaged (right-hand side) spectral estimates. (a), (b) FFT (rectangular window): (c). (d) AR order 16 (fpe); (e), i f ) MLM order 25.
C. Spectrum Estimation 0.
.5
1.
1.5
Time (sec)
2.
0.
.5
1. 1.5 Time (sec)
2.
(C) (d) Fig. 7 . Frequency-time curves for a shear layer (a) Burg AR peak tracker. ( b ) MLM peak tracker. ( c ) FFT peak tracker, id) FFT first moment.
phase locked loop systems [ 171. For such broad-band signals, a tracking algorithm can be used which takes into account the past history of the velocity signal and disallows sudden changes in the Doppler frequency. However, such algorithms are sensitive in that one looses the track easily when the velocity changes rapidly, a situation encountered physiologically. The MLM provides the most stable peak tracker. Although it is associated with a small underestimation of the fluid velocity, it is nearly as stable as the first moment estimators. Moreover, AR peak trackers have been found to be very sensitive to the choice of model order, even in the range 3-7, whereas the MLM tracker was more robust. Thus, the MLM peak tracker is preferable if the nature of signal (broad or narrow band) is not known.
It is desirable to display the spectrum of the signal because clinicians take a qualitative impression of spectral broadening into account during diagnosis [ 181. Under any experimental conditions, but especially in a shear layer, inspection of single spectra reveals large fluctuations in the spectral shape with time, even though the flow is steady. Spectral averaging greatly reduces these fluctuations and yields spectra which are consistent with the velocity profile independently of the spectral estimation method. However, this technique is difficult to implement reliably in pulsatile flow. The question thus arises as to which algorithm yields the most stable spectral estimates. Fig. 8 shows averaged and single spectra for data recorded in the strong shear layer. The spectra on the right [Fig. 8(b), (d), and (f)] are ensembles over 64 records. The figures on the left [Fig. 8(a), (c), and (e)] are 12 individual spectra plotted over one another. The three ensembled spectra have similar shapes. The unensembled realizations of the FFT spectrum in Fig. 8(a) show the widest variation. The AR spectra in Fig. 8(c) show less variability, and the MLM spectra in Fig. 8(e) show the least variability.
D A V I D er
< I / , :
595
MODERN SPECTRAL ANALYSIS TECHNIQUES
The MCM behave similarly to the Burg algorithm for both single and averaged spectra, independently of the model order, and hence the MCM spectra are not displayed. Furthermore, low orders for the modem methods result in unrealistic spectral shapes. Windowing the FFT does not significantly change either single or averaged spectra. For the AR methods, all orders above the FPE yield results similar to the one plotted. The FPE indicated orders between 12 and 18. For the MLM, an order of 25 was ultimately chosen because it yielded spectra consistent with the velocity distribution inside the sample volume for every degree of shear encountered. Higher orders were statistically less stable and computationally less efficient. IV. DISCUSSION The AR and MLM techniques have been shown to be superior to Fourier methods under nearly all experimental conditions when a peak tracking algorithm is used to find the Doppler frequency, provided the model order is chosen appropriately. An order of 5 performed well under all flow, noise, and signal bandwidth conditions investigated. Applying a Hanning window to the signal did not improve the behavior of the FFT. The performance of all the peak trackers was dependent on the bandwidth of the signal, i.e., the range of velocities present in the sample volume. The MLM yielded the most robust peak tracker in that it was the least sensitive to changes in model order and signal bandwidth. The Burg algorithm appears to be more robust than the MCM for peak tracking with the Doppler signals. For model orders above 15, numerical ill conditioning sometimes occurred for the MCM, and the MCM would also result in peaks which did not represent Doppler frequencies. These peaks were sometimes stronger than the actual Doppler shift peaks and led to spikes in the velocity-time curves from the peak detector. When the first moment technique was used to compute the velocity, the autoregressive and Fourier methods were virtually identical. The computation was nearly independent of the model order for the AR algorithms, and again, nonrectangular window shapes did not improve the FFT results. The first moment approach was barely sensitive to the signal bandwidth, but strongly affected by the presence of noise. The relative performance of the peak tracking and first moment methods depends on the level of noise in the signal and on the range of velocities present in the sample volume. The peak trackers were much more accurate (both in terms of bias and variance) than the first moments when noise was present in the signal. This was observed for both colored and white noise. However, they became unstable when the spectrum of the signal was broad, unlike the first moments. For spectral estimation, the MLM algorithm appeared to be more stable than the other approaches. Its use should be particularly rewarding in the analysis of pulsatile flow
signals, where the estimation of the velocity distribution inside the sample volume relies on single spectral estimates because reliable spectral averaging is difficult. It must be remembered, however, that the MLM spectrum is not a true power spectrum. The area under the spectrum between two frequencies is not the power in that frequency band. However, the shape of the spectrum still provides insight into the velocity field within the sample volume. Care must be taken if these spectra are to be used to extract specific power related indices. Comparison of the results presented here to those of Vaitkus and Cobbold [9] must be made cautiously because different model orders and different algorithms were used in the computations. Also, the spectra here apply to pulsed Doppler ultrasound whereas those of Vaitkus and Cobbold model continuous wave spectra. However, the lower spectral variance for the MLM method over both Fourier and AR methods is notable. The difference in optimal AR model order for the two studies, 16 for this study, 8 for Vaitkus and Cobbold reflects the broader nature of the spectra studied here. Questions still exist in regard to the applicability of the modern spectral analysis techniques to more complicated spectra, such as might be obtained from a stenosed artery. The following strategy for processing pulsed Doppler ultrasound signals can be proposed from the above observations: display simultaneously the MLM order 5 peak tracker, a first moment estimate, and a higher order MLM spectrum (an order of 25 seemed to be satisfactory for all the flow conditions studied). This procedure allows an estimate of the bandwidth of the signal, and thus a choice of the best suited velocity estimator (the peak tracker should always be preferred except for very broad spectra). This strategy is being proposed with currently available Doppler instruments equipped with FFT processors [ 191. It is clear that such a strategy can be optimized by using modem spectral analysis algorithms instead of FFT’s, since the AR yields first moments identical to those of the FFT, while the MLM peak tracker is always superior to the FFT tracker when the order is kept below IO. Since the MLM spectrum is computed from the AR parameters, it is not unrealistic to compute the first moment of the Burg spectrum while displaying the MLM spectrum. V. CONCLUSION Two important parameters influencing the accuracy of velocity estimation from pulsed Doppler ultrasound signals as obtained with four spectral analysis algorithms have been studied. It was shown that the optimal algorithm is dependent on the flow (bandwidth) and the signalto-noise conditions. Robustness considerations lead to the choice of a combined display of MLM order 5 peak tracker and high order (25) MLM spectrum. An AR first moment estimator can also be displayed if broad-band signals are expected. The results may be clinically significant since accurate and robust velocity estimators are needed in order to de-
596
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 38. NO. 6. J U N E 1991
rive reliable diagnostic ~ l visual~ inspection ~ , of spectra as well as quantitative studies on spectral shape (e.g., spectral broadening index) require stable and accurate spectral estimators.
REFERENCES W . K . George and J . J . Luniley. "The laser Doppler velocimeter and its application to the measurement of turbulence," J . Fl~tidM c d i . . vol. 60, pp. 321-362. 1973. S. M. Kay and S. L. Marple, "Spectrum analysis--a modern perspective." Proc. /EEL-. vol. 69. pp. 1380-1419. 1981. L. J . D'Luna and V . L. Newhouse. "Vortex characterization and
1171 D. P. Giddens and A . M. A . Khalifa, "Turbulence nieasureinents u'ith pulsed Doppler ultrasound employing a frequency tracking method.'' Ultrctsoir~ltlM o d . Biol.. vol. 8. pp. 427-437. 1982. 1181 W . J . Zwiebel. "Spectrum analysis in carotid sonography." Ultrci,sourid M r d . Biol.. vol. 13. pp. 625-636. 1987. 1191 W . F . Voyles. S. A . Altobclli, D. C. Fisher, and E. R. Greene. " A comparison o l digital and analog inethods of Doppler spectral analysis for quantifying flow." Ulrro.sorirrtl M e d . Biol.. vol. I I pp. 727733. 10x5. I
Jean-Yves David received the engineering degree from the University of Technology of Comidentification by ultrasound Doppler," Ulrrct.soiiic. /i~~ei
