Application of autoregressive methods to multigate spectral analysis

Ultrasound in Med. & Biol., Vol. 26, No. 4, pp. 585–592, 2000 Copyright © 2000 World Federation for Ultrasound in Medicine & Biology Printed in the USA. All rights reserved 0301-5629/00/$–see front matter

PII S0301-5629(00)00145-9

● Original Contribution APPLICATION OF AUTOREGRESSIVE METHODS TO MULTIGATE SPECTRAL ANALYSIS GABRIELE GUIDI, LEONARDO CORTI and PIERO TORTOLI Electronic Engineering Department, University of Florence, Via S. Marta, 3, 50139 Florence, Italy (Received 16 September 1999; in final form 6 January 2000)

Abstract—Multigate analysis is known to be capable of detecting accurate blood velocity profiles from human vessels. Experimental systems so far presented in the literature use time-domain frequency estimations and, more recently, the fast Fourier transform (FFT) for real-time analysis of Doppler signals from multiple range cells. This experimental study is aimed at evaluating the application of an autoregressive (AR) method (Burg algorithm) to multigate Doppler analysis. Both in vitro and in vivo results were collected with a commercial Duplex scanner coupled with a prototype multigate unit developed in our laboratory. The same multigate signals are, thus, processed according to both the FFT and the Burg algorithms. The related spectral and maximum frequency profiles are reported and statistically compared. AR, implemented with the Burg algorithm, is demonstrated to be a way to perform multigate spectral analysis with reduced spectral variance, suitable for maximum velocity profile extraction through a simple threshold. © 2000 World Federation for Ultrasound in Medicine & Biology. Key Words: Ultrasound Doppler, Velocity profiles, Flow profiles, Multigate, Specrtral analysis, FFT, Autoregressive spectral estimation, Burg’s method, Spectral variance.

extracted, rather than the single mean value attainable from correlation methods. In particular, it may be significant in extracting the peak frequency because this may be reliably converted to velocity, after ISB compensation (Tortoli et al. 1995; Winkler and Wu 1995). In Doppler systems, spectral estimation was traditionally achieved with the fast Fourier transform (FFT). Although faster and computationally robust in comparison with other algorithms, FFT presents drawbacks like spectral leakage due to the windowing that is inherent in finite length data sequences and, when ensemble averaging cannot be used (Bartlett 1948; Welch 1967), large spectral variance. For this reason, parametric methods have been proposed as possible alternatives to the FFT. In these methods, the signal is modeled as the output of a linear system characterized by a rational transfer function, whose input is a stationary random process. This approach provides a better frequency resolution and reduces spectral leakage, avoiding the use of windowing (Proakis and Manolakis 1992) at the expense of a major computational power and the need of imposing some assumptions on the signal behavior. In this study, the feasibility of introducing a parametric method for flow profile estimation is discussed. After reviewing the main features of autoregressive (AR)

INTRODUCTION The possibility of investigating blood-flow dynamics through the detection of velocity profiles intercepted along an entire ultrasound (US) beam has been demonstrated by several authors (Bassini et al. 1979; Hoeks et al. 1982; Robinson et al. 1994; Tortoli et al. 1985). Such profiles are generated by multigate systems that, in most cases, provide a single estimate for each group of Doppler samples related to a given depth. Auto or crosscorrelation algorithms are typically employed to estimate the mean velocity within the related range cell (Bonnefous and Pesque 1986; Embree and O’Brien 1985; Kasai et al. 1985). As an alternative, Doppler signals can be analyzed by estimating their whole spectrum. The Doppler spectrum, although its width can be influenced by factors not dependent on the flow dynamics—such as intrinsic spectral broadening (ISB)— can be approximately seen as a histogram of velocities, where the amplitude of each frequency component is related to the number of red cells producing that frequency. From the instantaneous spectrum, different spectral indexes can be Address correspondence to: Gabriele Guidi, Electronic Engineering Department, University of Florence, Via S. Marta, 3, 50139 Florence, Italy. E-mail: [email protected] 585

586

Ultrasound in Medicine and Biology

methods, this article describes an experimental procedure allowing the application of both AR and FFT methods to the same sets of multigate Doppler data. Examples of the application of this analysis to in vivo data are presented. Finally, in vitro results show that a higher accuracy in the definition of velocity profiles may be obtained with the parametric approach. AUTOREGRESSIVE SPECTRAL ESTIMATION A fundamental problem in the estimation of a signal power spectral density resides in the practical constraint of analyzing it over a finite time interval. If the signal is statistically stationary, the estimation can be affordably accomplished by taking a data sequence long enough to make time windowing effects neglectable. On the other hand, Doppler signals originating from blood flow are typically nonstationary, and the data sequence length for generating a spectral estimate is determined by the rapidity of the time variations in the signal. AR spectral estimation has been demonstrated to be suitable for short data records (Proakis and Manolakis 1992) and, for this reason, it has been successfully applied to Doppler spectral estimation (Vaitkus et al. 1988; Ahn and Park 1991), although spectral estimation accuracy may be greatly reduced in case of low signal-to– noise ratio (Proakis and Manolakis 1992). AR methods are, in general, based on modeling the signal sequence x(n) as the output of a linear system characterized by a frequency transfer function of the form: 1

H共 f 兲 ⫽

冘a e

(1)

p

1⫹

k

⫺j2␲fk

k⫽1

whose input is the sequence w(n). The power spectral density of x(n) can, therefore, be expressed as: P xx共 f 兲 ⫽ 兩H共 f 兲兩 2P ww共 f 兲,

(2)

where Pww(f) is the power density spectrum of the input sequence w(n). If this sequence is a zero mean white noise with variance ␴w2, then the power density spectrum of the input signal becomes:

P xx共 f 兲 ⫽ ␴ w2兩H共 f 兲兩 2 ⫽

冏

␴ w2

冘a e p

1⫹

k

k⫽1

冏

⫺j2␲fk

2

,

(3)

Volume 26, Number 4, 2000

and x(n) is called an AR process of order p (Proakis and Manolakis 1992). Spectral estimation through AR modeling consists in the determination of the ak parameters according to a proper algorithm. A linear system whose parameters are represented by the autocorrelation matrix of the signal must be solved. To evaluate the autocorrelation matrix, time samples of the signal before and after the time window considered have to be generated, and the various methods employ different approaches to extend the sample set. We selected the Burg method, as also suggested by other authors (Vaitkus et al. 1988; Ahn and Park 1991; Marasek and Nowicki 1994), because it offers good performance and, because of its recursive structure, is suitable for computer implementation and can be optimized for real-time applications. Moreover, the Burg algorithm offers good performance with respect to the problem of spectral leakage. Because our multigate system does not perform any high-pass filtering of input data, the output spectra are characterized by clutter components up to 60 dB stronger than the Doppler components. Burg’s method proposes an autocorrelation matrix estimation that does not use zero padding, as done, for example, by the Yule–Walker algorithm, but guesses the time samples to maximize the entropy of the process (Burg 1967). This approach, thus, reduces the possible spectral leakage of the strong clutter components, which might even cancel the weak informative signals. Whichever is the chosen algorithm, one of the problems encountered when dealing with rational spectral estimations is selecting a proper model order, p. Although an “optimum” model order could be selected through different techniques (Akaike 1969, 1974; Parzen 1974; Rissanen 1983), for our purposes, we selected a fixed model order p for every estimation. A good tradeoff was experimentally shown to lie in the range 10 –12 (Keeton et al. 1997; Schlindwein and Evans 1990). Kaluzynski (1989) suggested that, to model properly the Doppler spectrum in all phases of the cardiac cycle, the order should not be lower than 8. However, it has to be noticed that all the referenced papers consider spectral estimation starting from real signals. In our case, as described in the next section, because IQ demodulated signals are employed, we processed the complex extension of the Doppler signal. The AR spectral estimation is, therefore, evaluated starting from complex signals and this allows reduction of the model order to one half the order used for real signals because the complex conjugate poles pairs of the real analysis are not modeled in the complex approach (Kay 1978; Kay and Marple 1981). For this reason, in our experiments, we used p ⫽ 5 with complex signals, equiv-

AR methods and multigate Doppler ● G. GUIDI et al.

587

range of the A/D converters has to be large enough to record without clipping the small Doppler components originated by flow, superimposed to the huge low-frequency clutter signal due to slow movements of the vessel walls. In our equipment, a 12-bit device was chosen. The A/D sampling frequency can be set by the operator up to 8 MHz, equivalent to a minimum rangecells spacing of about 100 ␮m. Because no integration of the received samples is performed, the actual spatial resolution of the system (i.e., the sample volume dimensions) depends on the radiofrequency (RF) front-end (duration of transmission pulse, receiver bandwidth) as well as on the transducer focusing properties. In the following experiments, the Doppler equipment was set to provide sample volume lengths of 1–2 mm and sample volume widths lower than 2 mm at the range of interest. Fig. 1. FFT vs. Burg spectral estimation of the Doppler signal collected in vitro from a steady laminar flow with a 1-mm long range cell located near the tube axis.

alent to a model order of 10 with real signals, as suggested in the literature. MATERIALS AND METHODS Multigate system The multigate system used in our experimental test was described in detail in previous papers (Tortoli et al. 1996, 1997). It is based on a processing unit that can be connected to any Doppler equipment providing as external signals the pulse-repetition frequency (PRF) and the outputs of the phase-quadrature demodulators (I and Q components). The PRF is used to regenerate the sampling burst needed to perform a proper selection of 64 range-cells along the insonating beam. The analog-to– digital (A/D) section converts the I/Q signal components not subjected to any high-pass filtering. The dynamic

Experimental procedure The multigate system capability of performing the FFT processing in real-time is exploited to set the system and to position the US probe in the most appropriate way. We can, in fact, observe the so-called spectral profile on the real-time display and move the probe laterally until the wider extension in range is obtained, indicating that the beam axis crosses the vessel axis. When possible (e.g., in vitro), the probe-to–tube distance is changed until the probe focus is made coincident with the vessel axis. After the probe position has been fixed, the optimal range gate position and range cell spacing are set to cover an adequately large number of space bins. The PRF is finally selected to provide the needed extension in frequency free of aliasing artefacts. When all these parameters have been adjusted, a command from the computer keyboard starts the acquisition of the raw time-domain data into the computer RAM. The data can be saved on the hard disk, organized in 256 frames each composed of 64 rows of 64 I and 64

Fig. 2. 3-D instantaneous spectral profiles of a steady laminar flow generated in vitro, calculated with spectral estimation based on (a) FFT and (b) AR Burg.

588

Ultrasound in Medicine and Biology

Volume 26, Number 4, 2000

Fig. 3. Instantaneous spectral profiles obtained in vivo by coding the estimated power spectral densities to grey levels. CCA in diastolic phase: (a) FFT; (b) AR Burg.

Q samples. Finally, they are postprocessed through a set of MATLAB procedures for evaluating spectral profiles based on FFT and AR estimation techniques. In vitro flow phantom Some experiments were produced with an in vitro flow system using a solution of distilled water and 50 ␮m Sephadex particles. The concentration used was 1 g of Sephadex per L of water, equivalent to 0.1% in weight. This solution flows from an upper to a lower reservoir through a 1-meter long straight tube, with internal diameter of 10 mm and 1 mm of wall thickness. Such a tube stays immersed in water along all its length, so that the US probe can be positioned far enough from the tube inlet to avoid flow profile distortions. As specified in the literature (Caro et al. 1978), such distance must be at least 40 diameters, equivalent to 40 cm in our case. To have a known velocity profile, laminar flow conditions were set by imposing a Reynolds number lower than 2000 (Bird et al. 1960). With the given

diameter, and considering the cinematic viscosity of water, such limitation on the Reynolds number involves a maximum axial velocity of 40 cm/s. In our in vitro experiments, it was set at about 18 cm/s. RESULTS Figure 1 shows the instantaneous spectra estimated by analyzing the same Doppler signal with the FFT and AR Burg algorithms, respectively. The signal was backscattered from a range cell located near the tube axis, in steady laminar flow conditions. The figure highlights the expected smoothing effect associated with the AR Burg spectral estimate (Vaitkus et al. 1988). When all the 64 range cells are simultaneously measured, spectral profiles like those shown in Fig. 2 are obtained. The higher smoothing introduced by the AR Burg spectral estimation is clearly visible in these 2-D plots. Examples of in vivo spectral profiles originated by the same temporal data are shown in Figs. 3 and 4. In

Fig. 4. Instantaneous spectral profile obtained from analysis of CCA during systolic peak. (a) FFT based spectral estimates; (b) AR Burg.

AR methods and multigate Doppler ● G. GUIDI et al.

Fig. 5. Effects of clutter spectral leakage. FFT instantaneous spectral profiles evaluated by the same data producing Figs. 3 and 4 when rectangular windowing is superimposed: (a) Diastolic phase corresponding to Fig. 3; (b) systolic peak corresponding to Fig. 4. The white vertical trace present in both figures is due to a software zeroing of the first frequency bin, made to avoid an automatic grey scale too dark for allowing a proper profile display.

Fig. 6. Ensemble average of 256 instantaneous spectral profiles obtained from in vitro laminar flow with (a) FFT; (b) AR Burg algorithms; and (c) the superposition of the FFT and AR Burg spectra related to the central depth (28 mm) show that, at a threshold level of about ⫺30 dB, equivalent peak frequencies are estimated.

589

590

Ultrasound in Medicine and Biology

Volume 26, Number 4, 2000

Fig. 7. Statistical properties of the peak profiles obtained on the same set of data originating (Fig. 6) at a threshold of ⫺25 dB. Both spectral estimation methods were employed on instantaneous spectral frames: (a) FFT; (b) AR. In each figure, the standard deviation bars are plotted together with the peak frequency profile obtained from the average spectral frame of Figs. 6a and b, respectively.

these cases, related to the analysis of the common carotid artery (CCA) of a young male volunteer, the power spectral densities are coded to grey levels, as in conventional spectrograms. In Fig. 3, obtained during diastole, narrow spectra are involved because of the low blood velocity yielding to quasilaminar flow conditions. All the spectral power results, therefore, concentrated in few frequency bins and the Doppler spectrum emerges well above noise independently of the chosen algorithm. Figure 4, reporting the spectral profiles evaluated at the systolic peak, shows more irregularities on the FFT spectra (Fig. 4a) than on the AR spectra (Fig. 4b). The difference between this behavior and the previous one, may be attributed to the spreading of the same spectral energy over a larger number of frequency bins. It is worth observing that, while the AR-based spectral profiles were obtained from the raw Doppler data, the latter were weighted with a Blackman–Harris window to produce the FFT-based spectral profile. This window is necessary to reduce the spectral leakage effect due to the presence of a huge clutter in the signals to be processed. If windowing is omitted, FFT-based spectral profiles like those in Fig. 5 are obtained. The leakage effect is less evident in the diastolic phase when the arterial walls are approximately motionless (Fig. 5a), whereas the useful signal is almost destroyed by clutter at the systolic peak (Fig. 5b), when the vessel movement is stronger. It may be observed that, although the ⫺6 dB sample volume length was here set at about 2 mm, the clutter is so large that its contributions corresponding to the “tails” of the strong echo from the first wall are visible also in the center of the vessel. To get a quantitative estimation of the performance

obtainable with FFT and AR Burg methods applied to multigate analysis, a large amount of data was collected through in vitro experiments related to steady laminar flow conditions. A reference spectral frame was first obtained by averaging 256 subsequent FFT-based spectral profiles (Fig. 6a). Averaging was also applied on 256 sets of AR Burg parameters, producing a group of reference parameters that was used to evaluate a second reference spectral profile (Fig. 6b). It may be observed that the reference profiles obtained in the two cases, although not identical, are nearly equivalent. A reference “peak” profile was then extracted by fixing a threshold level 30 dB below the maximum power of the FFT-based spectral profile because, at such a level, the two peak profiles extracted from the FFT and AR spectra turn out to be the closest to each other. This is highlighted by Fig. 6c, showing the AR and FFT ensemble averaged spectral estimates, originating from a range cell located in the center of the vessel. The upper frequency cross-point between the two spectra lies at approximately ⫺30 dB from the maximum power spectral density. The 256 FFT- and 256 AR-based instantaneous peak frequency profiles were evaluated at the ⫺30-dB

Table 1. Average bias and standard deviation over the 64 range cells analyzed in Fig. 6 with a ⫺30 dB threshold Bias

FFT AR Burg

Standard deviation

Absolute (Hz)

%

Absolute (Hz)

%

⫺25.6 ⫺31.8

⫺1.5 ⫺1.9

67.2 42.1

3.9 2.5

AR methods and multigate Doppler ● G. GUIDI et al.

Table 2. Average bias and standard deviation over the 64 range cells analyzed in Fig. 6 with a ⫺15-dB threshold Bias

FFT AR Burg

Standard deviation

Absolute (Hz)

%

Absolute (Hz)

%

⫺43.8 ⫺60.1

⫺2.7 ⫺3.7

71.1 47.0

4.4 2.9

threshold level specified above, and compared with the FFT-based peak profile, assumed as the experimental “gold standard.” Standard deviation bars, computed for each of the 64 range cells, are plotted in Fig. 7a and b for the two processing techniques. The 64 values obtained in this way for bias and standard deviation were then averaged to provide the two statistical indexes reported in Table 1. These indexes, which are related to the profile as a whole, indicate slightly better performance of the Burg algorithm in terms of standard deviation, and an opposite situation when bias is considered. These features are maintained when a different threshold level is considered, as shown in Table 2, obtained by fixing the threshold level at 15 dB below the peak power of the reference spectral profile. DISCUSSION AND CONCLUSION Autoregressive methods have been so far considered as a possible alternative to the FFT in the analysis of a single Doppler signal. Their use was mainly suggested by the need to overcome the FFT limitations in terms of accuracy and precision of the spectral estimate (Vaitkus et al. 1988). As a consequence of the rational form of an AR spectrum, a spectral smoothing is also achieved. This paper has shown that, by applying the Burg algorithm in multigate analysis, such smoothing involves not only each depth at subsequent time slots, but also different depths at a same time. A regularization in space at a specific time is, thus, obtained. This feature can become of practical interest when single parameters like the peak frequency have to be extracted from a full spectrum. Simple thresholding methods (where “peak” is the frequency at which a fixed threshold is overcome) can be used because spectrum smoothness involves less variance (more precision) of the related parameter. On the other hand, because of the regularization of each instantaneous spectral profile, the corresponding frame turns out to be low-pass filtered, and possible details may be lost. When the full spectral profile is considered, no spectral analysis method appears definitely preferable to the other because both FFT and AR allow the flow

591

behavior to be clearly identified. An ideal system should allow both processing methods to be used. The processing load associated with Burg’s algorithm turns out to be finally compatible with the use of appropriate architectures based on digital signal processors (Guidi et al. 1998), especially after considering the limited model order needed for analysis of complex Doppler signals. Acknowledgements—This work was supported by the EC Grant BMH4-98-3782 (DOLPHINS Project). The authors thank Francesco Guidi for the precious technical assistance and Dr. Federico Frosali for some theoretical discussion about the various autoregressive spectral estimation techniques.

REFERENCES Ahn YB, Park SB. Estimation of mean frequency and variance of ultrasonic Doppler signals by using second-order autoregressive model. IEEE Trans Ultrason Ferroelec Freq Cont 1991;38(3):172– 182. Akaike H. Power spectrum estimation through autoregression model fitting. Ann Inst Stat Math 1969;21:407– 419. Akaike H. A new look at the statistical model identification. IEEE Trans Automatic Contr 1974;19:716 –723. Bartlett MS. Smoothing periodograms from time series with continuous spectra. Nature 1948;161:686 – 687. Bassini M, Dotti M, Gatti E, Pizzolati P, Svelto V. An ultrasonic non-invasive flowmeter based on cross-correlation techniques. Ultrasonics Int Proc 1979;NA:273–284. Bird RB, Stewart WE, Lightfoot EN. Transport phenomena. New York: John Wiley & Sons, 1960. Bonnefous O, Pesaque P. Time-domain formulation of pulse-Doppler ultrasound and blood velocity estimation by cross-correlation. Ultrason Imaging 19865;8:75– 85. Burg JP. Maximum entropy spectral analysis. Proceedings of the 37th Meeting of the Society of Exploration Geophysicists, Oklahoma City, 1967. In: DG Childer, ed. Modern spectrum analysis. New York: IEEE Press, 1967. Caro CG, Pedley TJ, Scroter RC, Seed WA. The mechanics of circulation. Oxford, UK: Oxford University Press, 1978. Embree PM, O’Brien WD Jr. The accurate ultrasonic measurement of volume flow of blood by time-domain correlation. IEEE Ultrason Symp Proc 1985;2:963–966. Guidi F, Guidi G, Ricci S, Atzeni C, Tortoli P. High speed parallel processing of biomedical ultrasound signals. Proceedings of the Second European DSP Education and Research Conference, Paris 1998:327–330. Hoeks APG, Reneman RS, Peronneau PA. A multi-gated pulsed Doppler system with serial data processing. IEEE Trans Sonics Ultrason 1982;28:242–247. Kaluzynski K. Order selection in Doppler blood flow signal spectral analysis using autoregressive modelling. Med Biol Eng Comput 1989;27:89 –92. Kasai C, Namekawa K, Koyano A, Omoto R. Real-time two-dimensional blood flow imaging using an autocorrelation technique. IEEE Trans Sonics Ultrason 1985;32:458 – 464. Kay SM. Maximum entropy spectral estimation using the analytical signal. IEEE Trans ASSP 1978;26:467– 469. Kay SM, Marple SL Jr. Spectrum analysis—A modern perspective. Proc IEEE 1981;69(11):1380 –1419. Keeton PIJ, Schlindwein FS, Evans DH. A study of the spectral broadening of simulated Doppler signals using FFT and AR modelling. Ultrasound Med Biol 1997;23:1033–1045. Marasek K, Nowicki A. Comparison of the performance of three maximum Doppler frequency estimators coupled with different spectral estimation methods. Ultrasound Med Biol 1994;20(7): 629 – 638. Parzen E. Some recent advances in time series modeling. IEEE Trans Automatic Contr 1974;19:723–730.

592

Ultrasound in Medicine and Biology

Proakis JG, Manolakis DG. Digital signal processing. Principles algorithms and applications. 2nd ed. New York: Macmillan Publishing Company, 1992. Rissanen J. A universal prior for the integers and estimation by minimum description length. Ann Stat 1983;11:417– 431. Robinson TM, Lee JW, Roberts VC. Detection of early atherosclerosis by analysis of ultrasonic Doppler signals produced by mural flow disturbances. Med Biol Eng Computing 1994;32-6:702–703. Schlindwein FS, Evans DH. Selection of the order of autoregressive models for spectral analysis of Doppler ultrasound signals. Ultrasound Med Biol 1990;1:81–91. Tortoli P, Guidi F, Guidi G, Atzeni C. Spectral velocity profiles for detailed ultrasound flow analysis. IEEE Trans Ultrason Ferroelec Freq Cont 1996;43:654 – 659. Tortoli P, Guidi F, Guidi G, Berti P, Atzeni C. A novel instrument for high-resolution ultrasound flow-imaging. Int J Imaging Syst Technol 1997;8:438 – 443. Tortoli P, Guidi G, Newhouse VL. Improved blood velocity estimation

Volume 26, Number 4, 2000 using the maximum Doppler frequency. Ultrasound Med Biol 1995;21:527–532. Tortoli P, Manes G, Atzeni C. Velocity profile reconstruction using ultrafast spectral analysis of Doppler ultrasound. IEEE Trans Sonics Ultrason 1985;32:555–561. Vaitkus PJ, Cobbold RSC, Jhonston KW. A comparative study and assessment of Doppler ultrasound spectral estimation techniques part II: methods and results. Ultrasound Med Biol 1988;14:661– 672. Welch PD. The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short modified periodograms. IEEE Trans Audio Electroacoust 1967;AU-15: 70 –73. Winkler AJ, Wu J. Correction of intrinsic spectral broadening errors in Doppler peak velocity measurements made with phased sector and linear array transducers. Ultrasound Med Biol 1995;21: 1029 –1035.