Determination of tissue motion velocity by correlation interpolation of pulsed ultrasonic echo signals

ULTRASONIC

IMAGING

12,

84-98

(1990)

DETERMINATION OF TISSUE MOTION VELOCITY BY CORRELATION INTERPOLATION OF PULSED ULTRASONIC ECHO SIGNALS P.G.M. de Jong,’ T. Arts,l~s A.P.G. Hocks,’

and R.S. Reneman*

Departments of Biophysics1 and Physiology* Cardiovascular Institute Maastricht University of Limburg P.O. Box 616 6200 MD MAASTRICHT, The Netherlands

Correlation interpolation is introduced as a method to determine the displacement of moving biological tissue on the basis of a sequence of ultrasonic echo signals. The echo signal is sampled along the echo depth with approximately 4 samples per average high frequency period. Sampling in time occurs with the pulse repetition frequency. The necessary information is extracted from a crosscorrelation function between successive signals, which is modelled using four parameters. The parameters are estimated from five calculated correlation sums and the shift with maximum correlation is determined. In contrast to existing techniques, the performance of this method is determined mainly by the number of samples used, while the ratio of the number of samples in depth and time is irrelevant. Using 64 samples at a signal-to-noise power ratio of 10, the standard deviation of the error in the determination of the shift in depth is 0.08 sampling intervals. As in many other methods, the width of the aliasing interval equals the mean frequency period. 0 1990Academic Press, Inc. Key words:

Correlation; Doppler; digital sampling; tissue; ultrasound; velocity.

echo; heart; signal analysis;

1. INTRODUCTION In the clinical environment cardiac contractile performance is usually evaluated by means of two-dimensional ultrasonic echocardiography (2DE). This technique provides a sequence of cardiac images at a rate in the order of 25 images per second. Since the method is based on detection of ultrasound reflections, the images show reflecting surfaces and scattering reflectors most clearly. As a result, cyclic changing cardiac geometry is deliniated by the moving epi- and endocardial con-

sAuthor 0161-7346/90 Copyright All rights

to whom correspondence

should be addressed.

$3.00 0 1990 by Academic Press, of reproduction in any form

Inc. reserved.

84

TISSUE VELOCITY

BY CORRELATION

INTERPOLATION

tours. The pattern of changes in cardiac geometry during the cardiac cycle provides information about cardiac performance [1,2]. For instance, in case of coronary artery disease, part of the myocardium performs poorly. The area and the degree of function loss may be estimated by means of echocardiography by analysis of the time course of wall thickness in the region of interest. Such analysis is complicated and subject to significant artifacts. One of the major problems in evaluating local cardiac function is associated with the difficulty of assessing motion of tissue in parallel with the reflecting wall surfaces. Most myocardial fibers are directed parallel to the wall. Hence, fiber shortening should be assessed from gradients in tissue velocity along a direction in which velocity measurement is difficult and inaccurate. Applying image analysis techniques on 2DE images, time dependent myocardial deformation of the cardiac wall can be studied by following reflecting texture images of the wall cross section [3-71. These methods function reasonably well in following global motion of the heart, but assessment of small differences in motion over a short distance is hampered by the limitations in spatial resolution of the 2DE image. The poor spatial resolution can be attributed to the use of envelope detection of the high frequency echo signal as used in the formation of 2DE images, thereby discarding phase information. Spatial resolution in motion detection can be improved by processing the high frequency ultrasound signal directly. Tissue velocity and, hence, tissue displacement can be determined with ultrasonic Doppler methods used regularly for the determination of blood flow velocity. Existing techniques employ the phase detector [B-l 01, the instantaneous frequency detector [l 1 ,121 or the complex autocorrelation detector [13,14], applied to the quadrature Doppler signal. The phase detector and the instantaneous frequency detector require a simple algorithm which can be performed in real time during data acquisition [15]. All these techniques operate on the principle that the echo signal contains a dominant frequency, shifting in phase with respect to a stable reference frequency. Between successive echo signals the phase shift at a certain depth is determined from a limited number of samples per signal. The number of samples that can be used meaningfully is limited by the finite correlation length of the high frequency ultrasound signal. When increasing the bandwidth of the signal, this correlation length shortens. Besides the use of pulsed Doppler methods, correlation techniques [16-201 can be used to determine tissue motion. In this approach the shift between two successive signals is determined by crosscorrelation of the signals. In contrast with the Doppler techniques, detection of the shift between successive signals is not limited by a possibly large bandwith, because larger bandwidth does not influence the correlation length between successive echo signals. In fact, the shift can be determined more accurately if more samples per signal are used. As a result, the resolution in depth is low whereas the resolution in time is high. The primary disadvantage of existing correlation techniques is the need of quite high computing effort. Generally, calculation of the crosscorrelation function requires a time consuming convolution. In case of relatively long signals, this convolution may be performed most efficiently by using forward and backward Fast Fourier transforms in succession [22]. In this article a method is presented which combines several advantages of the Doppler techniques and the correlation techniques. It is based on finding the signal shift with maximum crosscorrelation by the use of an interpolation algorithm. The analyzed high frequency echo signal may have a wide frequency band, and the advantageous noise suppression properties of crosscorrelation may be used. As in analyzing Doppler signals, the amount of calculations is limited because a

DE JONG ET AL

few evaluations of the correlation function in the vicinity of the expected maximum are sufficient to estimate the position of the correlation maximum. Furthermore, the number of samples used in echo depth and time can be choosen freely and independently. Of course, the more samples are used, the less the influence of noise is to be expected.

2. CORRELATION

INTERPOLATION

METHOD

2.1 General theory In the present study, the correlation interpolation method is introduced as a method to determine the relative shift between successive signals by analyzing a number of successive signals simultaneously. While periodically transmitting an ultrasound pulse, the echo signal received represents reflection intensity as a function of depth. When the reflecting structure moves with velocity vounder an angle awith the ultrasound beam, the successive echo signals shift gradually in echo depth (Fig. 1). The received signal y depends on two variables being depth x and time t: y(x,t)

= s(x-vt) + n(x,t)

with v=2

(1)

vocos a

The function s(x-vt) represents the reflectivity as a function of depth, which is moving in time with velocity v. The term n(x,t) represents noise which is uncorrelated with the function s(x-vt). A two-dimensional autocorrelation function is defined by: Ryy(X V = ~5, t Iyfx, t) y&+x, t+T)J

(2)

The operator E,r indicates expected value by random selection of x and tin the xt-area of interest. The variables X and T denote shift in depth and time, respectively. Assuming absence of correlation between the functions s(x-vt) and n(x,t), combination of Eqs. (1) and (2) results in: RyyfX

with and

?J

=

Rss(X

T) + &dX,

Rss(X, T) = Ex,t {s(x-vt) Rnn(X

7)

=

Ex,t

(3)

T)

{n(x,t)

s(x-vt+X-VT)] n&+X

t+Vl

Under the condition of ergodicity of the signal s(x-vt), the expectation on the basis of random selection of two variables x and t may be replaced by the random selection of one variable u=x-vt. A signal is defined to be ergodic if the statistical properties of the signal do not change with the variable at the abscissa [21]. The values of the function y(x,t) in Eq. (2) are available in a grid of sampling points with sampling intervals Ax and At for depth and time, respectively. Thus, the

86

TISSUE VELOCITY

BY CORRELATION

Fig. 1

depth

INTERPOLATION

Schematic representation of received echo signals. f f denotes the pulse repetition frequency and f, the depth sampling frequency. Displacement of signals due to motion is indicated by the broken line.

autocorrelation function &,(X,T) can be estimated numerically for discrete values of Xand T, to be defined by integers land J, respectively: . RYY(xJ T, = imax ’ jmax ‘y j=l

. Ikay i=l

y(x+X,

t+T)

(4)

with: x =iAx and t =jAf X=lAx and T= JAt i, j, I, J are integers

The numbers imaxand jmaxdefine the number of samples used for averaging in depth and time, respectively. The total numer of samples used equals the product of imax and jmax. Eq. (4) is used to estimate the autocorrelation function. The larger the number of samples used, the better the estimate will be. In the correlation interpolation method in the vicinity of the correlation maximum, the function Ryy(X,T) is approximated by a model. This model is based on Eq. (3) and contains only a few parameters. By calculating the estimate of Ryy(X,T) for a

87

DE JONG ET AL

few points, the parameter obtained.

values are solved and an estimate of the parameter

According to signal analysis [21, 221 the autocorrelation power density spectrum Y(w,T) of the function y&t) by:

function

R&X, T) = F -’ M-o> T) I

v is

is related to the

(5)

where the operator F indicates Fourier transformation, and o denotes frequency in echo depth. Assume the power density spectrum S(w) of s(u) is bandlimited and rectangular with center frequency oC, bandwidth wb, and power density So. Then Eq. (5) yields for R&X,T):

Rs,(X, T) = R,,(U) = F -' {S(o) I = T so cos(6&gJ)

sin (wbU/zL mbu/2

with U = X-VT

Furthermore, it is assumed that the power density spectrum of the noise term n(x,t) with fixed value of t has the same shape as the power density spectrum of the signal component s(x-vr). The similarity in shape is likely to be introduced by the bandpass filter characteristics of the ultrasound transducer used, which filter both signal and noise. As a result of the similarity of the signal and the noise spectrum, the noise term in Eq. (3) may be written as:

2 Kx(XT)

for T=O

RnnfXT,

=

0

(7)

for TAO

where No is the power density of the noise component. For T f 0, the noise component of the autocorrelation function equals zero because at different values of t the related signals are separately bandpass filtered, which implies absence of mutual correlation in the noise components of these signals. Using Eqs. (6) and (7), the autocorrelation function R&X,T) in Eq. (3) is found to be:

(1 + ‘$J Rss(X,T) R,#,T)

for

T = 0

=

(8) RssfX V

with Rss(X,T) autocorrelation

for T f 0

being defined according to Eq. (6). This model description of the function requires five parameter values,& No, wc, Wb, and v. The

88

TISSUE VELOCITY

BY CORRELATION

latter values may be obtained by calculating and solving the related set of equations.

2.2 Signals with narrow If the signal s(x-vt) plified to:

MXV

frequency

R&X,T)

INTERPOLATION

in five points using

Eq. (4)

band

has a narrow

frequency

band (ob
= : so cos(0, U)

(9)

Using Eq. (4), the autocorrelation function R,,(X,T) is estimated in 5 points being (O,O), (Ax,O), (O,At), (+Ax,At). According to Eqs. (8) and (9) Ryy(X,T) is written in terms of the parameters Ob, o,, SO, No and v. It holds:

Ryy(QO)

=T

(So + No)

(1W

RyyfAx 0)

= T

(So + No)cos(oc Ax)

(1Ob)

RyyfO~Af~

= T Socos(-co, v At)

(1w

Ryy@W)

= T so cos(w, Ax - oc v At)

(1W

Ryy(-Ax,At)

= 2 so cos (-wc Ax - wc v At)

(1W

This set of five equations contains four solvable unknowns. In order to solve parameter v, we have used the following calculation scheme. Using Eqs. (10a) and (1Ob), it follows:

oc Ax = arccos[Ryy(Ax,O)/R,,(O, 0)]

(11)

Subtracting Eq. (10e) from (10d) results in: T SOsin(oc Ax) sin (wc v At) = [Ryy(Ax,At) - Ryfi-Ax,At)]/2

(12)

The factor sin(oCAx) is calculated using the result of Eq. (11). Using Eqs. (12) and (1Oc), it is found:

oc vAt = arctan2([Ryy(Ax,At) - Ryy(-Ax,At)]/2, sin(wC Ax) Ryy(O,Af,))

89

(13)

DE JONG ET AL

The function arctan2(y,x) delivers the polar angular coordinate associated with the Cartesian x and y coordinates. Finally, the parameter v is found by division of the results of Eqs. (13) and (11) and multiplication by Ax/At. This calculation procedure is quite safe with regard to division by zero and invalid function arguments. The term R,#,O) denotes a sum of squares and is always positive as long as a signal is present. The argument of the arccos function in Eq. (11) is always between -1 and +l. The sampling interval Ax should be choosen so that sin(oC Ax) is in a range well beyond zero. The latter term is maximal if the center frequency is sampled approximately 4 times per period along the x-direction. To avoid aliasing in the arctan operation of Eq. (13), the angle (CI.I~vdt) should be within the range of -7~to n radians. Finally, according to the definition of v in Eq. (1) the velocity component in the direction of the transducer is found by halving the value of v.

2.3 Signals with wide frequency

band

In case of a larger bandwidth of the signal s(x-vt), the solution for v given by Eqs. (9) through (13) is not exact anymore because of a limited validity of Eq. (9). Nevertheless, the calculations can still be performed as if the signal had a narrow bandwidth. In figure 2 for signals with a flat, bandlimited frequency band, the relation between the exact solution on the basis of Eq. (8) and the estimate on the basis of the calculation scheme, as described by Eqs. (9) through (13), is shown for ratios of absolute bandwidth to center frequency of 0, 1 and 2, respectively. Since the central frequency is assumed to contain 4 samples per period, the highest frequency is sampled properly in all examples. From figure 2, it may be concluded that even for a relatively large bandwidth, being equal to the center frequency, the velocity estimator performs quite well. Because real echo signals may have quite a

RBW=O

RBW=

1

RBW=2 ---

Fig. 2

-3

-2

-1

0

1

2

displacement [sample points]

3

Theoretical analysis of the performance of the correlation interpolation algorithm when processing rectangularly bandlimited signals with a ratio of absolute bandwidth to center frequency (Relative Bandwidth) of 0, 1 and 2, respectively. The larger the bandwidth, the greater the discrepancy between calculated and real displacement will be.

TISSUE VELOCITY

BY CORRELATION

INTERPOLATION

large bandwidth, it is preferable the signals are preprocessed with a bandpass filter with a bandwidth less than or equal to the center frequency. By doing so, errors in the estimation of the velocity due to nonlinear characteristics are small (Fig. 2).

3. SETUP

OF

3.1 Bandlimited

SIMULATIONS

test signals

Echo signals as obtained by ultrasonic echographic equipment are simulated by bandlimited noise. To this purpose, noise ho(x) with a Gaussian amplitude distribution is generated with a random number generator, and bandpass filtered with a 6 point convolution bandpass filter (1 ,l,-2,-2,l ,l). This filter has a central frequency of 0.22 times the sample frequency (yielding signals with approximately 4.5 samples/period) and a bandwidth to center frequency ratio of 1.O. For testing purposes, successive signals hi(x) are made by shifting the noise signal ho(x) linearly in time. The number s represents the shift between successive signals, as expressed in number of sampling intervals. The shifted signal hi(x) is obtained by a combination of integer shifting and interpolation of ho(x) before bandpass filtering using: hi(x) = COS (3 hg(x+pAx) + sin (yj hg(x+(p+l)Ax) with p = in@)

and

(14)

r= js-p

By this interpolation method, the power of the signal is independent of the shift. Signals generated on the basis of different random number sequences are not correlated. Using this property, noise is obtained in precisely the same way as the signal. A noise signal of this kind is added to the signal with signal-to-noise power ratios of 10 and 1, respectively. Different scans have different noise components while the signal components are shifted versions of the same basic signal.

3.2 Comparison with other methods To compare the correlation interpolation method with other methods used for velocity measurement, we implemented two Doppler techniques. These techniques determine the frequency shift of the ultrasonic echo signal, due to displacement of the tissue, with respect to the transmission frequency. The Doppler frequency shift fdis related to the transmitted frequency f, and the tissue velocity vg as: f = 2vo cos d C

a fc

(15)

where cdenotes velocity of sound, and a is the angle between the motion direction and the ultrasound beam. We implemented the instantaneous frequency detector

91

DE JONG ET AL

[11,14] and the complex autocorrelation detector [13,14]. They both use the in-phase x&j) and quadrature yd(j) signal components, as obtained by demodulation of the j-th received high frequency ultrasound signal with the transmission frequency. The demodulation is performed using imax samples in depth. The instantaneous frequency detector estimates the Doppler frequency shift by averaging the following value over jmax signals in time:

$0) = &[

arctans

- arctan’%]

(16)

The Doppler frequency shift $0) is expressed as a fraction of the pulse repetition frequency. Alternatively, the complex autocorrelation detector estimates the Doppler frequency shift by:

In Eq. (17), averaging over jmax signals in time has been performed before performing the arctan operation. In both detectors, the arctan operation is implemented with the arctan function as used in Eq.(13). For narrow band signals, the frequency range of both detectors equals -0.5 to 0.5 times the pulse repetition frequency. All three methods are applied to the same set of test signals. In order to investigate the sensitivity to noise, the shift between signals is determined at signal-to-noise power ratios of 10 and 1, respectively. In order to evaluate the behavior of the different methods with respect to shape of the sampling area, the number of sampling points has been fixed to 64. So, the following sizes have been used: 2x32, 4x16, 8x8, 16x4 and 32x2 samples in depth and time, respectively. The complex autocorrelation detector and the instantaneous frequency detector cannot use less than 4 samples in depth per signal, because demodulation has to be performed over at least one signal period of the reference signal. This limitation does not hold for the correlation interpolation method. In order to evaluate the sensitivity to noise when changing the number of sampling points in depth and time, the error introduced has been assessed for the Doppler methods as well as the correlation interpolation method. For this evaluation the sampling areas used are the 9 combinations of 4, 8 and 16 samples in depth with 2, 8 and 32 samples in time. Generally, the Doppler methods are applied to signals with a relatively narrow bandwidth. The center frequency in depth is locked to the demodulation frequency and is generally a quarter of the sampling frequency. Thus, the calculated displacement between successive echo signals is the product of Doppler frequency shift and sampling frequency, divided by the center frequency of the signal. The correlation interpolation technique delivers the depth shift between successive echo signals directly. The echo signal may have a large bandwidth and the center frequency does not need to be locked to the sampling frequency.

92

TISSUE VELOCITY

BY CORRELATION

INTERPOLATION

4. RESULTS In figure 3, for signals with a bandwidth equal to the center frequency, the displacement as calculated with the correlation interpolation algorithm is shown as function of the real displacement. In absence of noise, the relation is quite linear up to the displacement where aliasing occurs. This aliasing occurs around half the periodicity of the center frequency of the signal. For displacements close to the aliasing point, the method seems less precise. The effect of the large bandwidth on the linearity is little, especially for small displacements. When noise is added, the location of the aliasing points, the slope and the linearity appear to be unchanged. Under the extreme condition that noise power equals signal power (Fig. 3, right panel), scattering of the measurements becomes large, but the method still works reasonably well. In table I, the sensitivity to noise is compared for the instantaneous frequency detector, the complex autocorrelation detector and the correlation interpolation method. Linearity and aliasing properties of all three methods are similar. However, the correlation interpolation technique appears to be less sensitive to noise. In this respect, the complex autocorrelation detector performs better than the instantaneous frequency detector. When increasing the number of echo signals used per measurement, all methods show a similar improvement in performance. The standard deviation does not vary significantly differently from the reciprocal of the square root of the number of samples used. In figure 4, the three methods are compared with respect to variation of the shape of the sampling area, while the number of samples used remains constant at 64. The results are obtained for sampling areas of 4x16, 8x8, 16x4 and 32x2 samples in depth and time, respectively. If depth is limited, the complex

a

,,,., -3-2-i

b

,,/,, 0

12

C

,,I

3

-3-z-i

real

0

displacement

12

[sample

,. 3

./, -3-2-l

I

,.

0

1

2

points]

Fig. 3 Calculated versus real displacement of bandlimited signals, as processed with the correlation interpolation method. An averaging sequence of 32 points in depth over 2 successive signals was used. (a) absence of noise; (b), (c) signal-to-noise power ratio of 10 and 1, respectively.

93

3

DE JONG ET AL

Table I. Standard deviation of the estimation of the displacement when processing bandlimited signals with three different methods at various signal to noise ratios. The standard deviation is calculated within the interval of -1.4 to 1.4 sample intervals displacement, because this interval appears to be free of aliasing in all cases investigated. Instantaneous Frequency Detector SNR=lO SNR=l

xt area

Complex Autocorrelation Detector SNR=lO SNR=l

Correlation Interpolation Algorithm SNR=lo SNR=l

4x2 4x8 4x32

0.33 0.25 0.18

0.77 0.37 0.49

0.35 0.20 0.07

0.63 0.37 0.31

0.23 0.11 0.05

0.73 0.37 0.19

8x2 8x8 8x32

0.51 0.20 0.15

0.84 0.42 0.54

0.50 0.17 0.09

0.90 0.67 0.26

0.19 0.08 0.04

0.62 0.26 0.13

16x2 16x8 16x32

0.34 0.20 0.15

0.86 0.49 0.51

0.39 0.18 0.13

0.83 0.50 0.24

0.15 0.05 0.04

0.47 0.17 0.12

1

autocorrelation detector and the correlation interpolation algorithm perform similarly. The aliasing area of the instantaneous frequency detector is smoothened due to averaging many signals, while the slope of the linear part is somewhat smaller than unity. While the total number of. sampling points remains the same at greater depth and shorter time intervals, performance of the correlation interpolation method appears to be nearly unchanged. Under these conditions, the performance of the complex autocorrelation detector deterioriates considerably. However, the instantaneous frequency detector improves with respect to unity slope and aliasing effect, and finally performs similarly as the complex autocorrelation detector. In this respect, the correlation interpolation method is superior in motion detection when averaging over many samples in depth.

5. DISCUSSION In order to determine the velocity of a depth signal moving in time, the correlation interpolation method is evaluated using test signals obtained by band limitation of noise. The method needs 5 correlation sums to be calculated over the region of interest. The shape of the correlation function in depth and time is estimated. Maximum correlation is determined on the basis of this estimate. Comparing the performance with existing methods based on Doppler frequency shift analysis [14], the correlation interpolation method appears to be less sensitive to noise in the signal. Added noise does not cause a systematic error in the estimate.

94

TISSUE VELOCITY

4x16

-3-2-l

0

1

2

3

-3-2-l

BY CORRELATION

INTERPOLATION

8x8

16x4

0

1

2

3

real displacement

32x2

-3-2-l

0

1

2

3

[sample points]

Fig. 4 Results of processing bandlimited signals with the 3 methods using averaging areas of different shape. Results of (a) instantaneous frequency detector, (b) complex autocorrelation detector and (c) correlation interpolation algorithm versus real displacement, respectively. Areas of 4x16, 8x8, 16x4 and 32x2 samples in depth and time were used.

Although the algorithm is based on sinusoidal signals, the displacement can be determined accurately when processing bandlimited signals with a bandwidth up to the central frequency (Fig. 3). The calculated displacement is linear with the real displacement, implying that nonlinearities, as induced by a larger bandwidth, are moderate, as is also shown by theoretical analysis (Fig. 2). If the bandwidth of the signal in depth is larger than the mean frequency, the bandwidth should be reduced by appropriate bandpass filtering before further analysis. However, this is not likely to occur because most echo systems have a relative bandwidth between 0.5 and 1 [23] due to the filtering effect of the transducer used. For many applications, such a band limitation is not a severe disadvantage. Because the frequency band is quite wide with respect to the mean frequency, the band limitation may be performed by convolution with a relatively short impulse response consisting of 5 or 6 points. Due to the limitation in bandwidth, the signal has a periodic appearance, causing aliasing. By this phenomenon, a second

95

DE JONG ET AL

correlation maximum occurs at a distance of the order of the periodic length corresponding with the center frequency of the bandpass filter. Thus, the area, in which the correct correlation maximum is found, is limited to the range from -l/2 to l/2 times this periodic length. Linearity and aliasing properties of the instantaneous frequency detector and the complex autocorrelation detector are similar to those of the correlation interpolation method. However, when using many signals over a short depth range, the slope of the instantaneous frequency detector tends to decrease (Fig. 4). This underestimation, which has been reported before 112,141, is due to erroneous contributions of instantaneous frequencies to the result, especially for wide band Doppler signals. Smoothening around the aliasing point of the instantaneous frequency detector is due to averaging (Eq. (16)). The latter two disadvantages are not present in the performance of the complex autocorrelation detector. The instantaneous frequency detector as well as the complex autocorrelation detector are sensitive to a change of the shape of the area used (Fig. 4; table I). When a short depth interval is used over many signals, the performance of the complex autocorrelation detector is similar to that of the correlation interpolation method. However, when using a large depth range, the performance deteriorates, because the required good correlation with the demodulation frequency deteriorates with increasing signal length. After all, wide band signals have short coherence lengths, and when the depth interval used for summation exceeds this length, the quality of correlation diminishes. In practice, pulsed Doppler techniques are used over approximately 16 depth samples, which is only possible when the bandwidth of the signal is small. With the correlation interpolation method, the depth interval used is not limited by such a correlation length (Fig. 4; table I). The correlation interpolation method is based on estimating expected values of correlation products. According to the theory, the selection of the area in depth and time, where the individual products are obtained, is irrelevant as long as the signals may be considered ergodic [21]. This irrelevancy is confirmed by the finding that changing the shape of the area of summation in depth and time hardly influences the performance of the correlation interpolation method (Fig. 4). So spatial resolution and time resolution can be exchanged freely. Using the instantaneous frequency detector, analysis of N sampling events requires 0.25 N moderately accurate arctangent operations. These operations can be performed quickly, enabling real-time handling of the signals received [15]. The complex autocorrelation technique is more complex. The number of multiplications needed is N, and therefore the method needs more computational effort. The correlation interpolation method is even more complex and requires 5N multiplications. Nevertheless, computational effort is still linear with the number of sampling points used, and unlike techniques based on finding the maximum correlation, iteration techniques are not needed. Calculation of the whole correlation function can be performed most efficiently using fast Fourier transformation, which requires many more computational operations.

6. CONCLUSIONS

The correlation interpolation method estimates the displacement between two or more bandlimited signals in succession. The method is quite insensitive to noise in

96

TISSUE VELOCITY

BY CORRELATION

INTERPOLATION

the signal. The number of mathematical operations needed is linear with the number of sampling points used for analysis. The width of the interval free of aliasing is one period of the mean frequency. Within this interval the method appears to be linear and accurate. The method can simultaneously average the signal contributions in depth and time. Increase of the depth averaging interval and decrease of the time averaging interval, while the number of points used remains constant, does not significantly influence the performance of the algorithm. The correlation interpolation method performs better than the instantaneous frequency detector and the complex autocorrelation detector, especially when the number of depth samples used is large and noise is significant.

REFERENCES

[ll

Troy, B.L., Pombo, J., and Rackley, C.E., Measurement of left ventricular wall thickness and mass by echocardiography, Circulation 45, 602-611 (1972).

PI

Mastrototaro, J.J., Willcox, T.O., and Pilkington, T.C., The use of twodimensional echocardiograms in the detection of myocardial infarction in canines, IEEE Trans. Biomed. Eng. BME-32, 621-629 (1985).

[31 Mailloux, G.E., Bleau, A., Bertrand M., and Petitclerc, R.,Computer analysis of heart motion from two-dimensional echocardiograms, IEEE Trans. Biom. Eng., BME-34, 356-364 (1987). PI

Trahey, G.E., Hubbard, S.M., and von Ramm, O.T., Angle independent ultrasonic blood flow detection by frame-to-frame correlation of B-mode images, Ultrasonics 26, 271-276 (1988).

PI

Sai Prasad, R., and Srinivasan, T.M., An image processing method for cardiac motion analysis, IEEE Trans. Biom. Eng. BME-34, 244-247 (1987).

PI Dickinson,

correlation

R.J., and Hill, C.R., Measurement of soft tissue motion between A-scans, Ultrasound Med. Biol. 43,263-271 (1982).

using

VI Tristam, M., Barbosa, D.C., Cosgrove,

D.O., Nassiri, D.K., and Bamber, J.C., Hill, C.R., Ultrasonic study of in vivo kinetic characteristics of human tissues, Ultrasound Med. Biol. 12, 927-937 (1986).

PI

Barber, W.D., Eberhard, J.W., and Karr, S.G., A new time domain technique for velocity measurements using Doppler ultrasound, IEEE Trans. Biom. Eng., BME-32, 2 13-229 (1985).

PI

Grandchamp, P.-A., A Novel Pulsed Directional Doppler Velocimeter: The Phase Detection Profilometer, in: Ultrasonics in Medicine, pp. 137-l 43 (Excerpta Medica, Amsterdam, 1975).

[lOI

Nowicki, A., and Reid, J.M., An infinite gate pulse Doppler, Biol. 7, 41-50 (1981).

Ultrasound

Med.

[l l] Brandestini, M.A., and Forster, F.K., Blood Flow Imaging Using a Discrete-time Frequency Meter, in 1978 Ultrasonics Symposium Proc., pp. 348-352 (IEEE, New York, 1978).

97

DE JONG ET AL

[12] Hoeks, A.P.G., Peeters, H.H.P.M., frequency estimator for sampled BME-37, 212-220 (1984).

Ruissen, C.J., and Reneman, R.S., A novel Doppler signals, /EEE Trans. Biom. Eng.,

[13] Kasai, C., Namekawa, K., Koyano, A., and Omoto, R., Real-time dimensional blood flow imaging using an autocorrelation technique, Trans. Sonics Ultrason. W-32, 458-464 (1985).

two/EEE

[14] van Leeuwen, G.H., Hoeks, A.P.G., and Reneman, R.S., Simulation of realtime frequency estimators for pulsed Doppler systems, Ultrasonic imaging 8, 252-271 (1986). [15] Hoeks, A.P.G., On the Development of a Multi-gate Pulsed Doppler System with Serial Data-Processing, Thesis, Univ. of Limburg, The Netherlands (1982). [16] Bolon, P., and Lacoume, J., Speed measurement by cross correlation theoretical aspects and application in the paper industry, /EEE Trans. Acous. Speech Signal Processing ASSP-31, 1374-l 378 (1983). [17] Wear, K.A., and Popp, R.L., Theoretical characterization of myocardium contraction of ultrasonic echoes, IEEE Trans. Ultrason. 368-375 (1987).

analysis of a technique for the based upon temporal correlation Ferroelec. Freq. Cont. UFFC-34,

[18] Bonnefous, O., and Pesque, P., Time domain formulation of pulse-Doppler ultrasound and blood velocity estimation by cross correlation, Ultrasonic Imaging 8, 73-85 (1986). 1191 Knapp, C.H., and Carter, G.C., The generalized correlation method for estimation of time delay, /EEE Trans. Acoust., Speech Signal Processing, ASSP-24, 320-327 (1976). [20] lanniello, J.P., Time delay estimation via cross-correlation large estimation errors, /EEE Trans. Acoustics, Speech ASSP-30, 998-l 003 (1982). [21] Peebles Jr, P.Z., Probability, Random Principles (McGraw-Hill, New York,1 980). [22] Bracewell, R.N., The Fourier Transform York, 1978).

Variables,

and

and its Applications

[23] Hoeks, A.P.G., Ruissen, C.J., Hick, P., and Reneman, evaluate the sample volume of pulsed Doppler systems, Biol. 70, 427-434 (1984).

98

in the presence of Signal Processing, Random (McGraw-Hill,

Signal New

R.S., Methods to Ultrasound Med.