The analysis of blood velocity measurements by autoregressive modelling

Z theo~ BioL (1986) 120, 419-442

The Analysis of Blood Velocity Measurements by Autoregressive Modelling R. I. KITNEY AND H. TALHAMI

Engineering-in-Medicine Laboratory, Department of Electrical Engineering, Imperial College of Science and Technology, London SW7, U.K. AND

D. P. GIDDENS

School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332, U.S.A. (Received 19 June 1985, and in revised form 28 October 1985) Analysis of velocity disturbances arising in poststenotic flow is examined using the technique of autoregressive modeling. The essential elements of this method are described including relevant criteria for selecting the model order. Velocity data employed in the analysis are taken from in vivo measurements in the dog aorta, and the results indicate that the autoregressive method improves the resolution of coherent features in disturbed flow patterns. By applying homomorphic filtering to individual beats, the occurrence of organized structures convected from their origin in the shear layer is readily identified.

Introduction Blood velocity waveforms may exhibit significantly disturbed features in poststenotic flow, and within the past several years numerous noninvasive commercial Doppler ultrasound instruments which purport to analyze these disturbances in clinical situations have become quite popular in vascular laboratories. However, an understanding o f blood flow disturbances and their relationship to disease states is as yet far from complete. Indeed, poststenotic flow patterns are among the most difficult fluid dynamic phenomena to describe inasmuch as (i) they occur in an inherently pulsatile, as opposed to steady, flow; (ii) they are created at relatively moderate Reynolds numbers where asymptotic flow behaviour cannot be assumed, with transition during the cycle being the rule rather than the exception; and (iii) they are essentially nonstationary, waxing and waning during the cardiac cycle, not only in amplitude but also in frequency content. Khalifa & Giddens (1981) demonstrated the existence o f three types o f flow disturbances distal to a stenosis: a coherent structure arising just subsequent to the initiation of each cardiac cycle and hence termed a starting structure; a periodic oscillation arising in the poststenotic shear layer; and random disturbances which appear to possess typically turbulent energy spectra. Furthermore, depending upon the degree and shape of a stenosis and upon the Reynolds number, these structures may interact in such a manner as to render 419

0022-5193/86/120419+24 $03.00/0

© 1986 Academic Press Inc. (London) Ltd

420

R.I.

KITNEY

ET AL.

diflicult the separation of repeatable from random disturbances. Kitney & Giddens (1982) examined the methods of ensemble averaging and low pass filtering for obtaining the underlying periodic behavior in poststenotic velocity waveforms and concluded that, particularly when phase-shift averaging is employed to reduce phase jitter in the cycle initiation, ensemble averaging is by far the superior approach. Determination of the underlying waveform by low pass filtering produced decidedly inferior results. Having isolated the underlying waveform, the random velocity may be obtained by subtraction of the ensemble average from individual beats. Furthermore, although low pass filtering is inappropriate for securing the underlying waveform, it may be very useful to employ high pass filtering of individual beats to aid in the identification of coherent flow features which are higher in frequency than the fundamental heart rate. Consequently, the characterization of flow disturbances in a pulsatile flow may require several different approaches in an effort to analyze coherent and random features and to describe their nonstationary evolution. This paper concentrates upon analysis of coherent features in poststenotic flow, while a subsequent paper will consider the description of evolving random velocities. In particular, the use of autoregressive methods is studied in detail using a data base taken from hot film anemometer measurements of blood flow velocity distal to subtotal vascular constrictions in the descending thoracic aortas of mongrel dogs (Giddens et al., 1976). The varying degrees of stenoses studied, coupled with the use of the hot-film anemometer which offers excellent frequency response to velocity fluctuations (in contrast to Doppler ultrasound methods), provide a wide range of reliable poststenotic data with which to evaluate the usefulness of autoregressive methods for this appplication.

Basic Theory The classical approach to the description of system dynamics is by means of the Laplace transform or, if the frequency response is required, by the Fourier transform. This paper concerns itself with the application of a different type of spectral description based on linear estimation. The principles of these techniques will now be discussed in outline. Much of the original work on the development of autoregressive techniques was carried out by Box & Jenkins (1976). For this reason AR models are often referred to as Box-Jenkins models. Let us assume that a linear process can be described by the output of a linear filter whose input is white noise " a " . Hence, for the tth sample z,, define ~'t = at + a~raat-i + ' ~ 2 a t _ 2

+.

. .

or

~, = ~ %a,_j j=o

(1)

where ~, = zt - / z , g is the mean value of z,, and ~o = 1. Thus, z, is the deviation of

BLOOD

VELOCITY

421

MEASUREMENTS

the process from some datum, or from its mean if the process is stationary, and q,~, xt'2 etc. are weights. The white noise source has zero mean and a variance er2 and consists of a sequence of uncorrelated random variables. Since the variables at are uncorrelated, it follows that their autocovariance function has the property that ~cr 2, k = O "Yk= L 0 , k~O. Hence, the autocorrelation function Pk becomes I, rrk= 0,

k=0 k#0.

(2)

This result implies that a given output zt can be expressed as the sum of weighted previous outputs plus a forcing function at, that is, (3)

zt = a l 3,_ 1 + a2~.t-2 + . . . + at.

Now, equation (3) can be re-written for the special case of a finite set of weight parameters ~bl,..., ~b,. Letting :~, = y, equation (3) becomes y, = ~blY t - i + ~;b2Yt-2+... + ~ p y ,-p + a,.

(4)

The process defined by equation (4) is an AR process of order p. It follows directly that Yt = ~bty,-1 + at and y, = ckl t , - i + q52Yt-2+ at

are first and second order AR processes, respectively. Returning to equation (4) and multiplying through by Y , - k , we have Y t - k Y t = ~ b l Y t - k Y , - i + qb2Yt-kYt-2 + dP3Yt-kY,-3+. • . + c k p y t - k y t - p + y t - k a t

(5)

which gives, on taking expected values, yk=qblyk-1+~2Yk-2+...+dppyk_p

for k > 0

where y is the covariance. Dividing by Yo the autocorrelation form is then Pk = ~ ) l P k - -

1 + ~blPk-2+ ~ 3 P k - 3 + . . "~pPk--p"

(6)

Substituting k = 1, 2 . . . , p in equation (6) gives P l = (Jbl + (J~2Pl + .

• • + #~pPp-1

pp = qblPp_ 1 + ~ 2 P p _ 2 + . . . + ¢~p.

These equations are known as the Yule-Walker equations which, rewritten in matrix

422

R. I. K I T N E Y

ET AL.

form, are

LilL11 1111111 or pp = Ppqbp,where pp = the theoretical autocorrelations. Hence, 4~p= Pp-~Ppand it is clear from equation (6) that the coefficients of the original AR process are the solutions of the equation. The significance of this result from the viewpoint of spectral analysis lies in the relationship between the autocorrelation function of a time series and its frequency spectrum. Wiener established that the autocorrelation function of a signal and its power spectrum constitute a Fourier transform pair, and a similar result was established by Khinchin. Consequently, the relationship is often called the WienerKhinchin relationship. The Blackman & Tukey method (1959) provided a practical implementation of the relationship. A detailed survey of these techniques was given by Kay & Marple (1981). All these methods suffer from the fundamental disadvantage that the autocorrelation function is unknown for delays larger than the length of data under consideration. The assumption is that the autocorrelation estimates are zero outside the data interval. Ideally, the process being studied should therefore be stationary. Since a number of types of biological signals are "locally stationary", it is reasonable to estimate the autocorrelation function over a short period for which equation (6) comprises a limited number of terms. Under these circumstances the Yule-Walker equations and the Weiner-Khinchin relationship provide a reasonable basis for spectral estimation.

CALCULATION

OF THE SPECTRUM

OF THE AR MODEL

In order to understand how the PSD is calculated from the coefficients of the AR model obtained by the solution of equation (6), let us consider the basic system transfer function G(s) expressed as the ratio of two polynomials in s (where s is the Laplace transform operator). Hence three processes can be defined for the linear time-invariant case in terms of G(s): (a) the autoregressive (AR) process, represented by an all-pole model; (b) the moving average (MA) process, represented by an all-zero model; (c) the autoregressive moving average (ARMA) process, represented by a polezero model. In the AR process the input to the system, X, is taken as a white noise source of zero mean and known variance. The sampled data form of the system output can therefore be defined as

Y(z)=G(z)X(z)

(7)

BLOOD

VELOCITY

MEASUREMENTS

423

where 1

G(z) -

(8)

P

1 + ~. akz -k k~l

The frequency spectrum is derived from the definition of the z transform, so that Y(f) -

(9)

p

1+ ~

a k e -j2~fl~T

k=l

where IV(f) is the white noise source. THE

MAXIMUM

ENTROPY

METHOD

The methods of spectral analysis discussed thus far suffer from the difficulty that they all assume that the autocorrelation estimates outside the data are zero. In order to overcome this fundamental problem, a new approach called the maximum entropy method (MEM) was developed independently by Burg (1976) and Parzen (1968). The aim of the method is to obtain a spectral estimate of a time series based on a minimum number of assumptions. Box & Jenkins (1976) show that the AR family models a process as Gaussian. The entropy of a stationary Gaussian process (Edward & Fitelson, 1973) is given by H = 4j~N'

In S(f) d f

(10)

N

where fN is the Nyquist frequency and S(f) is the power spectrum subject to the constraints o(n) =

I?

s ( f ) exp ( - 2 ~ r f f n T ) dr,

n = - N , . . . ,N.

(11)

fN

Because of the minimum assumption requirement, when the first N autocorrelation lags are not precisely defined, equation (11) should be replaced by a constraint that bounds the distribution of values. Hence, maximizing the entropy of the process will provide a maximally noncommittal estimate of the power spectrum, with respect to unavailable information. The method developed by Burg corresponds to extending the known set of autocorrelation values beyond the maximum lag while maximizing the entropy at each stage. Returning to the time series previously discussed in equation (4), i.e. y ( 1 ) . . , y ( N ) with sample interval T. If these samples are assumed to arise from a pth order AR process, then ia

~,,=-

~., qb,,yn_,,+a(n) m=l

(12)

424

R.I.

ET

KITNEY

AL.

where as before, ~b,, are the AR coefficients and a(n) is a white noise source of variance o'~. Equation (12) therefore;constitutes a pth order linear prediction filter. The Burg algorithm uses the concept of forward and reverse filtering involving complex conjugation. The basis of this approach lies in Burg's statement (Burg, 1972) that the prediction error power can be calculated in this manner. If the prediction error " e " is defined as the ditterence between the predicted value of the nth sample y, and its actual value y(n), then for a pth order filter the forward prediction error can be written directly from equation (12) as P

ep=yn-~,=

E ~bp.,y,~-m

(13)

m=0

where 4~ = 1. Similarly, the backward linear prediction error b is defined by bp = ~ ~p.v~_p+m

(14)

m=O

for p < n < N - 1 . To obtain estimates of the AR (predictor) parameters, Burg minimized the sum of the forward and backward prediction error energies N=I

Nzl

E~-- E le~ol2+ E lbpf n=p

n~p

(15)

which leads to the result N--!

-2 (~)ii -

~

b i* -l,m-lei-l,m

(16)

lv-I

~ m=i

(Ib,_l.m-d2+le,-~.ml ~)

by setting dEp./ dqbp. = O.

Application to Velocity Data In terms of signal analysis the poststenotic velocity data under investigation are time series of relatively short duration; this is true irrespective of whether the analysis is carried out on individual beats or relationships between beats are being investigated. As a general rule Fourier spectral analysis requires a minimum of approximately five cycles of the fundamental frequency for reasonable resolution (Bergland, 1969). Furthermore, with Fourier methods the spectral content of a waveform can only be defined in terms of the power at harmonics of the fundamental frequency, so that if a true component lies between two harmonics, its power will be divided. Autoregressive spectra have a number of important potential advantages, the principal ones being that reasonable spectral estimates can be obtained with as little as one cycle of the fundamental frequency (Ulrych, 1972) and that the autoregressive spectrum is continuous. The improvements in spectral estimation which can be obtained from autoregressive methods are associated with the fact that a particular model will be more

BLOOD VELOCITY

MEASUREMENTS

425

parsimonious, provided certain criteria are properly met, than the equivalent Fourier description. The main criterion which is discussed in the literature is that of model order (Ulrych & Bishop, 1975; Kay, 1979; Kay & Marple, 1981). Ulrych & Bishop (1975) suggest that as a rule of thumb the model order can be taken as laying in the range 1 / 3 N < - p < - 1 / 2 N , and a similar criterion has been defined by Linkens (1979). An alternative approach is to use the Final Prediction Error (FPE) or Akaike Information Criterion (AIC) developed by Akaike (1969). A second important factor is sampling rate which was discussed in an earlier paper (Kitney & Giddens, 1982) and will now be considered in more detail. The velocity measurements employed as a data base for this study were obtained using hot film anemometry in the descending thoracic aorta of a dog as reported by Giddens et aL (1976), and the experimental techniques are described in that reference. The hot film sensor was mounted on the tip of a right-angle hypodermic needle and inserted into the aorta. Blood flow disturbances were induced by constricting the aorta proximal to the probe using a snare-type device whose diameter could be controlled. The cases studied ranged from 0% (no constriction) to 88% reduction in area, and several examples of the measured velocity waveforms are given in Fig. 1. The development of poststenotic velocity disturbances can be seen to increase as the degree of stenosis increases, with the deceleration phase being more sensitive to the creation of disturbances than the acceleration phase. (a} 1oo V

"G

(c) lOOr

..i.

.

,

l~,

. . . . .

I,I..

0L

FIG. 1. Velocity waveforms measured 2 cm distal to aortic constrictions for (a) 20, (b) 58 and (c) 88% area reductions.

The method of decomposing velocity waveforms into the underlying ensemble average and the disturbance velocity has been described previously (e.g. Khalifa & Giddens, 1978; Kitney & Giddens, 1982). Briefly, the R-wave of the ECG signal is employed as a fiducial point in sampling the velocity waveform. The ensemble average is then constructed from a desired number of beats, N. Although the ECG

426

R.I.

KITNEY

JET A L .

is reasonably reliable as a trigger for sampling the velocity, several factors contribute to an additional variability in beat-to-beat phase changes within the velocity waveform itself. The heart is not a perfectly repeatable mechanical pump, for example; and the onset, strength and duration of systole may each vary slightly. Respiration also affects the reproducibility of the waveforms and is a factor which is particularly evident in the hot film data employed in the present study. Consequently, in an effort to improve beat-to-beat phase alignment, a phase-shift averaging method is introduced to construct the ensemble average. The zeroth order average Uo(t) is computed directly from the sample data 1

N

Uo(t) = ~ ,~, ui(t)

(0-< t < P )

(17)

where u~(t) is the velocity waveform for the ith beat, P is the period selected for analysis and N is the number of beats forming the ensemble. Following this, each beat ui(t) is crosscorrelated with Uo(t) and a phase is determined for the shift which corresponds to maximum crosscorrelation. After each u~ beat is shifted accordingly, a new estimate of the ensemble average is computed and the entire process is repeated until no further improvements are obtained. We denote the final phaseshifted ensemble average by U(t). Next, the disturbance velocity is determined from

u:( t) = u,( t)

-

U( t)

(18)

for each beat. This disturbance, in principle, contains only the random velocity. In practice, however, this ideal cannot be achieved. For example, if high frequency shear layer oscillations occur during the deceleration phase, these will not always be initiated at precisely the same phase of the cycle even though the primary features of the waveform have been phase-shifted. Thus, the u'(t) disturbance velocity may contain some degree of coherent phenomena in addition to a purely random contribution. MODEL ORDER

Phase-shifted ensemble average velocity waveforms were obtained for all degrees of stenosis (Kitney & Giddens, 1983). The data for 88% occlusion are employed to demonstrate effects of model order on the spectra of the underlying waveform. Figure 2 illustrates the results obtained with 512 data points (8th order model) and the corresponding FFT amplitude spectrum. The sampling rate for the velocity waveform is 2560 Hz. Panel " a " shows the velocity waveform and panels "b" and "e" the spectra and z plane diagram for the AR model. It is clear from the figure (plot b, dotted characteristic), that the AR spectrum calculated from an 8th order model is a very inaccurate description of the signal. A similar conclusion is reached on the basis of the z plane pole pattern where only one of the 8 poles is close to the unit circle. As the model order is increased the situation is similar for 16 and 32 order models. A reasonable AR spectral estimate can only be obtained from high model orders which correspond to the number of points in the Fourier spectrum.

BLOOD VELOCITY

"5

427

MEASUREMENTS

(e)J

-

0-5

0

0'05 0"10 0"15 0'20

0

Time (sec)

1.0j~,,

E

-1-0 -0-5 l

20 40 60 80 100 Frequency (Hz) (c)

0-5 xi:0

Reel

FIG. 2. Phase shift averaged ensemble data (100 beats). (a) The velocity waveform. (b) 8th order autoregressive spectrum (AR) (dotted line). (c) z-plane diagram of the 8th order AR" model for 512 data points.

Figure 3 illustrates the Fourier and AR spectra for this case. While high order descriptions may be adequate for some purposes, the AR estimate in this case would be entirely inadequate if it is to be used parametrically, i.e. if a process, such as evaluating the degree of stenosis, is to be described by the AR coefficients.

.g°'-~g100

0

0

(o)

~

(b)

6000 4000

5000 2000

13

o

0

20

40

60

80 Frequency

0 (Hz)

20

40

60

80

FIG. 3. Power spectra of the 88% occlusion velocity data of Fig. 2a. (a) FTI" spectrum. (b) Burg AR spectrum--model order 256. Both spectra calculated from a 512 point data set.

REDUCTION

OF S A M P L I N G RATE

It might be inferred from these results that AR spectral estimates are apparently unsuitable for these velocity data. In the background theory section the representation of sampled data in terms of z transforms was described. Referring to the z plane the circumference of the unit circle is equivalent to the frequency axis in the s plane, hence the zero phase point is 0, to~, 2tos, 3to~, etc. and the 180 ° point is tos/2, 3toJ 2 , etc. The original sampling rate for the velocity data was 2560 Hz, so that the

428

R. I. K I T N E Y

ETAL.

180° point corresponds to 1280 Hz. A fundamental aspect ofautoregressive modelling is that the model will attempt to describe the entire frequency range (0-oJs/2). Therefore, if the sampling frequency is significantly in excess of twice the highest frequency of interest in the signal, the AR model will include an estimate of at best irrelevant signal information and, at worst, noise. It is therefore often important to reduce the original sampling rate by sub-sampling the original signal prior to modelling. Returning to the example of the 88% occlusion data, Fig. 4(a) illustrates the spectrum calculated from a 16th order model and 64 points, i.e. reduction by a factor of 8. The tos/2 point on the unit circle corresponds to 160 Hz; and as the maximum frequency in the ensemble averaged velocity data is 100 Hz, there is still adequate bandwidth. Comparing the AR spectrum with those calculated from the full signal bandwidth, it is clear that there has been a significant improvement in the spectral description, and the z poles are also sitting near to the unit circle. The effect of increasing the model order to 32 is illustrated in Fig. 4(b). It is clear from the figure that the AR spectrum now gives a better estimate of the data than that obtained by Fourier analysis.

(o) 1,C

1.

~g E

0"5

~'-1.0 x -0"5 x

20

40 60 80 Frequency(Hz)

t,O

-0"

-|.

100

Reol

(b)

10 0

ioI°' ~ 0

20

40 60 80 Frequency(Hz)

-

:>,

-1.0~

1"0

100 Reol

FIG. 4. A R spectra of the 88% occlusion velocity data of Fig. 2(a). In each case, the amplitude spectrum and z-plane diagram. (a) Model order 16, data decimated by a factor o f 8. (b) Model order 32. In both eases the upper trace is the F F T spectrum a n d the lower trace the A R spectrum.

These results show that although in theory a 256 order filter with a sampling rate of 1280 Hz should result in a good spectral estimate because the majority of the signal content (>90%) is below 100 Hz, only about 8% of the unit semicircle is devoted to velocity information. With a 160 Hz sampling rate, 62% of the semicircle

BLOOD VELOCITY MEASUREMENTS

429

contains the signal and 38% noise. When calculating AR spectral estimates it is therefore important to optimize the sampling rate as well as the model order. CRITERIA

FOR THE SELECTION

OF MODEL

ORDER

Although the choice of model order has been shown to be important, criteria for the selection of model order have not been discussed. This section of the paper will therefore consider the use of model criteria in the study of velocity data. The ensemble average velocity waveform consists mainly of coherent structures which are not removed during the averaging process. Hence, the ensemble average waveform is weakly stationary. Burg's Maximum Entropy Method (MEM) nevertheless can be applied successfully provided the right model order is chosen. For a weakly stationary time series of length N, low values of model order give insufficient resolution due to a smoothing effect of the true spectrum. Conversely, large values of p produce a peaky spectrum with statistically erroneous estimates (Haykin, 1979). We have studied the effect of the model order on the spectrum of the ensemble average velocity waveform and other simulated signals. Two model order criteria were used: the Final Prediction Error (FPE) and the First Zero Crossing of the autocorrelation function (FZC).

The final prediction effort ( FPE ) criterion The final prediction error is defined by Akaike (1969) as the expected variance of the prediction error when an autoregressive (AR) model fitted to a series of X ( n ) is applied to another independent realization of X(n). Hence, considering the time series X(n),n = 1, 2, . . . , N then FPE of X ( n ) = E ( X ( n ) - ~ ' ( n ) ) 2

(19)

where X(n) is the value estimated from the AR model. Akaike (1969) has shown that the FPE consists of two components: the first corresponding to the error of the best linear predictor for a given model order, and the second due to statistical deviation of 8p(m) from ap(m) where 8p(m) is the ruth estimated model coefficient for order p and ap(m) is the true coefficient. As the model order increases, the first term in the FPE decreases while the second increases for a given length of N of X(n). The result of Akaike's derivation can be stated as FPE (p) -

N+p+l

N-p-1

PM

(20)

where PM is the output error power of the filter. PM decreases with p while the other term which is due to the fluctuations of aM(m) increases with p. The FPE function has a minimum, Pmi,, for which, according to the criterion, p will be optimal, i.e. the optimal model order. This is only true if the stochastic process under observation is an AR process generated from a strictly stationary and mutually independent innovation. The relative FPE (RFPE) defined as (FPE),, (RFPE)p = (FPE)o

(21a)

430

R.I.

KITNEY

ETAL.

where N+I FPE0 - N - 1PM

(21b)

is the criterion ussed in this study. By way o f an example we applied the FPE criterion to a simulated 4th order stationary A R process. Figure 5(a) shows a realization o f the 4th order AR series

X ( n ) = alX(n - 1)-a2X(n - 2 ) - a3X(n - 3 ) - a 4 X ( n - 4 ) + e(n) where a 1 = - 0 . 2 3 7 3 ,

a2=+0-1940,

a3=-0-6312,

^

g

(22)

and e(n) is a

a4=+0-7293,

(b)

41 ~ 0 3 ~0"2 , "-E : ~o.1 o o

2-2 <{ 0,2

L

0,08

0 4 0.6 0 8 Time (secs)

0

(c)

10

20

30

50

40

Frequency (Hz) (d)

0.2

0-06 o 0.04

~. 0-02 0 0

Q

10

50

-1,1

I

20 30 4 Frequency (Hz)

50

i

0 0

I

50

I

i

30

40

50

Order 100 150 I

I

200

250

I

i

(f)

(e)

-1.!

20

Frequency(Hz)

Order 100 150 200 250 I

10

-0"5

-2,( -2,'

-1 C

~'-2.2 1.13

-1,5

~ -2.3 -2-4 -2.6

2

,

5

~

-2-0 -2-5

-2"7 FIG. 5. A realization o f a 4th order autoregressive model. Model order 256, from a 512 point data set. (a) The time series, (b) FFT power spectrum. (c) Burg-AR spectrum calculated from the time series. (d) AR spectrum calculated from the filter characteristics. (e) Final prediction error as a function of model order. (f) Mean prediction power as a function of model order.

BLOOD VELOCITY MEASUREMENTS

431

Gaussian bandlimited white noise source of variance o'o. 2 Comparing the power spectrum obtained using the FFT, Fig. 5(b), with the theoretical MEM spectrum, Fig. 5(d), calculated from the coefficients applying equation (16), it is clear that the FFT spectrum has some random fluctuations. The Burg MEM, Fig. 5(c), for p = 256, N = 512, overspecifies the model order resulting in a spurious spectrum with peaks not representative of the true frequency content o f the time series. The FPE (equation (21a)) gives a minimum at 4 (Fig. 5(e)) which is the correct model order for the process. Note that the FPE in this case is characterized by a sharp fall and then rise (Fig. 5(e)) and therefore is o f a form consistent with Akaike's criterion. The mean output power of the filter, Fig. 5(f), drops sharply for p < 4 and continues decreasing with higher model orders but at a slower rate. For p = 4 the spectrum is consistent with that predicted theoretically. The FPE for the 88% velocity ensemble average waveform of Fig. 2(a) reduced by a factor of 5, gives a minimum at p = 3. The corresponding spectrum is illustrated in Fig. 6. Note the overall smoothing of the spectrum and the large component at

v

4o

o.

(/3

0

20

40

60

80

Frequency(Hz) FIG. 6. AR spectrum of 88% occlusion velocity data of Fig. 2(a). Model order 3, using 512 data points.

zero frequency. It is therefore clear that the optimal order chosen by the FPE criterion has resulted in a bad estimate with very poor resolution. This result is typical for the application of the FPE criterion to ensembled averaged blood velocity waveform data. A significant improvement can frequently be obtained by manual selection o f model order around the FPE minimum. An alternative criterion based on the first zero crossing o f the autocorrelation function (the FZC criterion) gives a much improved estimate.

The first zero crossing criterion ( F Z C ) As shown in the theoretical section of the paper, the model coefficients are related to the autocorrelation coefficients by the Yule-Walker equations (equation (13)), the size of the autocorrelation matrix being equal to the model or filter order. It will be shown that the resolution of a harmonic in the MEM spectral estimate is a function of how well the extension of the ACF based on the Burg estimate matches the true ACF, which lies in the range + oo. As discussed previously, the accuracy of the prediction of the ACF depends on the model order; however, there is a minimum requirement for model order p.

432

R.I.

ET AL.

KITNEY

For a signal of n sample values the A C F can be calculated directly in the time domain over a m a x i m u m shift range of ± ( n - l) lags or + 1) seconds, where T is the sample interval. Given p A R coefficients the A C F can be extended recursively beyond this point. Hence, the A C F can be represented by

T(n -

IR=(n)

f°rlnl--
r=(n)=~_ k~=,ankr=(n_k)

for I n I > p

(23)

R,==

where true ACF, r,= = calculated ACF. For a sinusoid of frequency fo the A C F is a cosine with the same frequency and the expression in (23) is equivalent to multiplying the true A C F with a window which can be defined as w(n) = ~ 1.0

[ f(n)

for I.l<-p forlnl>p

(24)

where the shape o f f ( n ) is dependent on p and on the noise level in the original time series. A simulated single beat was used to investigate the effect of a window on the spectral estimate. The simulated beat is illustrated in Fig. 7(a) and consists of a "1 Hz sine wave to which Gaussian white n o i s e - - s t a n d a r d deviation 5 - - h a s been added. The section o f the beat used in the calculations is 0-45 of a cycle (sampling rate 64 Hz), the point o f primary interest being how effectively can the A R method resolve the fundamental frequency for this fraction o f a cycle. For the FPE minimum model order (p = 4), the Burg spectrum (Fig. 7(b)) has a primary peak at f = 0. The

601 40

~I0

~2

(~)

20 8

E

t\2o

4o

6o

so

~oo ~2o

-2oF

o 0"2

04 0.6 Time (secs) 1.0 0-5 t.L

<~

0'8

~

io

)

o -05 -1-0

FIG. 7. Half sine wave test signal with bandlimited random noise added. (a) The time series. (b) Burg AR spectrum--model order 128 256 points. (c) Autocorrelation function of waveform "a" together with 0-5 second unity window with exponential fit.

BLOOD VELOCITY MEASUREMENTS

433

true peak at 1 Hz is not resolved. The ACF of the signal (Fig. 7(c)) shows an exponential decay in the theoretical function as the number of lags is increased beyond 4. Referring to equation (23) this corresponds to the condition lnl > p where the observed ACF, r~,(n), can be described by the product of the true ACF, R~,, and a window function w(n). Here, r~,(r) = e -°~ cos toz

(25)

where e -°~ is the window function and ~" the time variable o f the ACF (Fig. 7(c)). In the frequency domain the window is convolved with the 1 Hz spectral component and produces resultant power in the zero frequency region (see Fig. 7(b)). As the model order is increased the unity range of the window expands (equation (24)), i.e. r~_~= R=, for n <- p, and the low frequency resolution improves concomitantly. Figure 8(a) shows the Burg spectral estimate for p = 16 from which it can be seen that the 1 Hz component has now been resolved. 60

(b)

I

(0)

4 0 i" Frequency (Hz)

2O

''2s ,0

~

ot

03 -20

-401

FIG. 8. AR spectrum of the test signal (Fig. 7(a)) decimated to 64 data points. (a) Model order 16. (b) Model order 32.

The basis of the first zero crossing criterion (FZC) is that p = 16 corresponds to the point where the ACF crosses the time axis, i.e. n = 16, and is also a function of sampling rate. Hence, the FZC criterion provides a reliable means of obtaining the minimum model order consistent with adequate low frequency resolution. Figure 8(a) shows the AR power spectrum calculated for a model order of 16. Referring to the figure the fundamental frequency is now clearly resolved in spite of the fact that the original time series is only 0-45 o f a cycle. With p = 32 (Fig. 8(b)) there is peak splitting due to a noise component close to f = 1 Hz, indicating that the model order is too large. The need to compromise between good low frequency resolution and modeling of noise components is of importance in AR spectral estimation. The previous example illustrates that increasing model order moves the " h a r m o n i c " poles towards the unit circle, increasing the unity segment of w(n). However, the noise poles will also move closer to the unit circle (Kay, 1979); and if the model order is too high, this will cause the appearance of "false" peaks, or peak splitting, in the spectrum. Hence, an optimum model order to produce good low frequency resolution should be high enough to achieve this objective while at the same time not be too high in order that problems with noise poles are avoided. This can be achieved by using the FZC criterion.

434

R. I. K I T N E Y

ET AL.

Figures 9(a)-(c) illustrate the use of the FZC in defining model order for velocity data from three different levels of occlusion 0, 40 and 88%, respectively. The original sampling rate was 2560 Hz but because the ensemble average velocity waveform comprises low frequencies only, the signal was decimated by a factor of 8, i.e. to a new sampling rate of 320 Hz. In all three cases it was clear from the original time series that the fundamental frequency was 5 Hz, and from the ACF's that the FZC occurred at n = 16. Consequently, applying the FZC criterion, the model order was set up to 16. Referring to Figs 9(a)-(c), there is a primary peak in all three spectra at 5 Hz. It should be noted that at the decimated sampling frequency the number of data points per cycle is 64; hence, as in the previous example, the FZC occurs after a quarter of a cycle.

40

~

L

o 20

co

-40

(b)i

(o)

oo _. frequency (Hz)

40

8.

60

requency (Hz)

20

40

60

(c)

} 4o "o 20 co

0 o

20

40 60 80 Frequency (Hz)

FIG. 9. A R spectra of ensemble average velocity waveforms, 512 data points decimated to 64 points. (a) 0% occlusion. (b) 40% occlusion. (c) 88% occlusion.

The application of the FZC criterion to a range of velocity data has shown that it is superior to the FPE criterion when good resolution of the underlying waveform is important, e.g. when studying the ensemble average velocity waveform. A second important point is that by following judicious application of decimentation and the FZC criterion, the AR power spectral estimates are consistently superior to the corresponding FFT estimates even when the duration of the time series is less than one cycle of the fundamental. It is important to note, however, that the accuracy of the AR estimate degenerates below one cycle and is heavily dependent on the portion of the cycle used (Ulrych & Clayton, 1976). A N A L Y S I S O F I N D I V I D U A L BEATS

The analytical examples considered thus far have focused on the ensemble average waveform for different levels of occlusion: equally important is the analysis of

BLOOD VELOCITY MEASUREMENTS

435

individual beats. The fundamental signal model is that an individual beat consists of an underlying component u(t) plus a higher frequency component, the disturbance velocity u'(t). As these two components are additive (equation (18)), the disturbance velocity o f an individual beat can be extracted from the overall waveform by subtracting the ensemble average waveform. An important limitation o f the approach is that both the simple and phase shift ensemble average techniques optimize the accurate representation of the low frequency content of the signal at the expense o f the high frequency content. Hence, when studying the high frequency content of individual beats, high pass filtering is likely to be a more appropriate alternative. A 100 Hz Kaiser filter was used for this purpose, and Fig. 10 shows the reduction of an individual beat (40% occlusion) into its low and high frequency constituents, c(t) and h(t), respectively. Considering the high frequency constituent h (t), the FFT power spectrum IH(to )l2, Fig. 10(d), illustrates that the frequency content of h(t) seems to comprise an

(a)l 5

5O

0

~

(b)

v

~ 25 I I

0 0

0-05

0-70

lO

0.15

0 0.05 TPme(secs)

0-~0

0~5

(c) 4C ~ |d

io

Frequency(Hz}

E -40

-lo L

0-05

.L

L

0-10 0.15 Time(secs)

FIG. 10. Analysisof 40% velocitydata, an individual beat. (a) The original time series. (b) Waveform "a" Iowpass filtered (Kaiser) to 100 Hz. (c) Residual time series of "a" followingfiltering, h(t). (d) FFT power spectrum of waveform "'c" h(t).

436

R.I. KiTNEY

ET AL.

underlying spectral trend with noise superimposed. In order to test this observation the Burg MEM spectrum was calculated according to the FPE criterion. As discussed earlier, the FPE minimum is designed to give the optimum model order for accurate description of the stationary components of a signal. Hence, from equation (10), the stationary component s(t), of the waveform h(t), can be defined in the z plane as

1

S(z) -

(26)

P

1 + ~ ak2-k k=l

where p is the model order defined by the minimum value of the FPE, and H ( f ) can be written as

H(f) = S(f)N(f)

(27)

on the assumption that the noise component approximates to a white noise source. Figure 11 (a) shows the Burg spectrum of the waveform of Fig. 10(c) calculated with a 20th order AR model, 20 being the FPE minimum. Comparing the FFT spectrum (Fig. 10(d)) with the Burg spectrum shows that the 20th order model is a good estimate of the underlying process. This is confirmed by the reconstruction of the AR model (Fig. 1l(b)) which is similar in structure to the original signal (Fig. 10).

I~,~

(o) I

Ib) 5

50o

~ooo

0 05

010

0.15

T~me(sec)

FIG. 11. (a) Burg AR spectrum of the waveform of Fig. 10(c). (b) Reconstruction of h(t) Fig. 10(c)) from the Burg realization.

In general, h(t) for individual beats can be modelled satisfactorily by this method if the model order is the minimum value of the FPE. Similar results were obtained for all the occlusion levels studied. Figures 12(a) and (b) are the original and reconstructed h(t) waveforms for a single representative beat at the 88% occlusion level. Although the Burg spectral estimate shows that h(t) comprises wideband noise and underlying spectral trend, the method is inappropriate for studying the process of interaction of these components in the time domain. Inspection of h(t) over

BLOOD

VELOCITY

(a)

10 A

t

5

E ~ 0 ~> - 5

437

MEASUREMENTS

(b) 10-

~o

-10 L 0-05

_L 0.10

I

0-05

0,15

Time (secs) FIG. 12. Original and reconstructed data, single beat.

h(t) w a v e f o r m s - - s e e caption

0-10 0.15 Time (secs)

of Fig. 11 for details--88% occlusion

many beats leads to the conclusion that around peak systole there is a region of peak activity, e.g. Fig. 10(c), where the general structure of the signal approximates to a type of DSB amplitude modulation with suppressed carrier. Thus, h(t) can be modelled as h(t) = s,(t), s2(t)

(28)

where s=(t) is equivalent to the carder signal and the signal s~(t)--the modulating signal--is a slowly varying waveform with a periodicity equal to the heart period. The waveform h(t) can be demodulated by homomorphic filtering (Oppenheim and Schafer, 1975). Referring to Fig. 13, the signal h(t) is first converted to log form,

x(n•) X

Log

~

Lineor system (filter)

FIG. 13. A block diagram o f the h o m o m o r p h i c filter procedure used to demodulate the waveform s h o w n in Fig. 14(a).

h(t)

then filtered in the frequency domain, low pass for sl(t), high pass for s~(t). The time series s~(t) and s2(t) are then obtained by retransforming and calculating their antilogs, i.e. exponential functions. Figure 14(a) illustrates an example of two beats of 40% occlusion velocity data prefiltered to remove components below 100 Hz. Figure 14(b) is the Fourier power spectrum of the function loge h(t). Signals 1 and 2, Fig. 14(c) and 14(d) are the demodulated signals s~(t) and s2(t).

438

ETAL.

R. I . K I T N E Y

(b)

(o) lO

"•0'4 o~

o

~" 0.2 o -10

# 0.2

0 0

0.4 0.6 Time (secs)

20

40 60 80 Frequency (Hz)

100

|

5 -

~

120

(d)

0.251-

,~ 2.5

o 0.2

0.4 0.6 Time (secs)

0

0.2

0.4 0.6 Time (secs)

FIG. 14. (a) Two beats of highpass filtered (100 Hz), 40% occlusionvelocitydata, h(t). (b) The power spectrum of h(t). (c) The modulating signal sl(t). (d) The carrier signal s2(t). From a fluid mechanics standpoint the importance of the method lies in its ability to clarify the structure o f the flow field at the measuring site, i.e. downstream of the stenosis. As a result of analysis over a range of velocity data a flow model has been developed which is based on the hypothesis that the major contribution to the disturbance velocity at higher levels of occlusion consists of organized structures o f discrete frequency which are convolved in time. The individual structures are taken to be of the form x ( t ) = A t 2 e -t/t° cos tot

(29)

where to is the wavefrequency and tc is a characteristic decay time. Figure 15(a) shows an example of a high pass filtered velocity for a single beat at 40% stenosis. Referring to the waveform, it can be seen that there are regions of relatively stable sinusoidal activity in terms of frequency, but with highly variable amplitude. If the signal is subjected to homomorphic filtering the low frequency envelope can be extracted (Fig. 15(b)). Referring to the figure, the points A, B, C, D and E define the amplitudes and times at which the organized structures pass the velocity probe. Hence, in this case the high pass filtered waveform of Fig. 15(a) contains six structures. The remaining high frequency content is shown in Fig. 15(c), and there is a period of apparently discrete frequency of oscillation (ca. 100 Hz) occurring between t = 0.10 and 0.15 sec. Figure 15(d) gives an impulse train approxi-

439

BLOOD VELOCITY MEASUREMENTS (a)

10

E

5

:=

0

3

4-

(b)

5

C

2

g 1

o

,-~J

0

0.05

Time (secs)

|

010 0-15 Time (secs)

0.20

(d)

~4

0.25

'5 0

> 0

-0,25 0.20

0

,t, , 1 1 , OO5

Time (secs) 100

0.10 0.15 Time (secs)

0.20

(e)

5O

~

0

0

~ -50 Time (secs)

FIG. 15. (a) A single, highpass filtered (100 Hz) beat taken from 40% occlusion velocity data. (b) Low frequency envelope of part (a). (c) High frequency content of part (a) after homomorphic filtering. (e) The result of convolving part (c) with a sequence of six impulses (d) of appropriate amplitude and occurrence time to form the reconstructed waveform.

mation for the amplitude function shown in Fig. 15(b). If this train and the frequency ¢ac = 100 Hz, derived from Fig. 15(c), are employed to the system described by equation (29), the reconstituted waveform shown in Fig. 15(e) results. While not identical to the original waveform, it can be seen that the main features are indeed preserved. Based upon previous experience with stenotic flows it is likely that the peak A in Fig. 15(b) is associated with residual disturbances, created during the previous cycle, being convected past the probe location as systole begins, while peak B arises from the starting structure described previously. The remaining peaks appear to be caused by the passage of large scale structures, such as shed vortices, created in the shear layer formed by the constriction. For example, if the Strouhal number (St = fd/Up where f is the frequency of vortex shedding, d is the constriction diameter and Up is the peak value of the ensemble averaged velocity) is calculated for the frequency associated with peaks C and D, the result is St= 0.57--a value in close

440

R. I. K I T N E Y

ETAL.

agreement with poststenotic Strouhal numbers reported by Cassanova & Giddens (1978) and for confined and free jet flows by Beavers & Wilson (1970). The high frequency content, represented by Fig. 15(c), is likely due to finer scale phenomena within the vortex structures. Conclusions

The present study demonstrates the need for an armoury of signal analysis techniques if maximum understanding of poststenotic fluid dynamic behaviour is to be attained. Traditional methods of time averaging and simple Fourier analysis are simply inadequate. Several important conclusions can be derived from the present study. (1) The use of phase-shift ensemble averaging is required to obtain the underlying velocity waveform in poststenotic flow. If low pass filtering is employed, either important cyclically-reproducible features are smoothed or random components are retained, depending upon whether the filter setting is "too high" or "too low". The result is a waveform which is not representative of the repeatable components of the underlying flow behaviour. Further with the variability in heart rate and phasic occurrence of flow patterns that is typical in biological systems, it is important to use the phase-shift averaging technique rather than a simple ensemble average. (2) Once the underlying waveform is obtained by phase-shift averaging, its frequency content is characterized better by autoregressive methods of spectral estimation than by Fourier methods. Although the waveforms are weakly nonstationary, the ability of AR methods to handle relatively short data records allows a reliable description of important frequency components. However, the questions of model order and data decimation must be carefully considered. The final prediction error criterion for model order selection often proves inadequate for blood velocity waveforms in that the model order selected by this method cannot identify the fundamental frequency present and it tends to smooth important higher frequency peaks. The first zero crossing criterion was shown to be preferable for this particular application. (3) When high frequency components of velocity disturbance are of interest, it was found that high pass filtering is a useful first step. This method allows both random and reproducible features to remain, but reduces the large amplitudes associated with the fundamental frequency, thus allowing a better focus upon phenomena of interest. The resulting high pass filtered signal is also nonstationary hence AR spectral estimation methods are likely to be preferable to FFT techniques when attempts are made to describe the spectral evolution during a cycle. (4) The nature of the HPF signal for moderate degrees of stenosis suggests that an important contribution to poststenotic flow disturbance arises from discrete frequency activity which is amplitude modulated. By employing homomorphic filtering to the HPF velocity it was shown that the nature of the 40% stenosis data is such that during a segment of the cycle the velocity can be modelled by a convoluting of individual structures of the form x ( t) = A t 2 e - ' / ' ,

cos tot.

BLOOD VELOCITY MEASUREMENTS

441

The frequency of arrival of these structures corresponds to that expected for vortex shedding from the constriction. Thus, the homomorphic filter approach allowed the recognition of an envelope of coherent activity immersed within the velocity disturbances. The origin of the discrete high frequency component is uncertain; however, we speculate that one possibility for its source is vibration of the vessel resulting from excitation by the poststenotic disturbances. If this is shown to be the case, it would be a very interesting coupling phenomenon between flow and the wall. Additional experiments in which wall motion itself is measured would be instrumental in this regard. In summary, the present study has yielded new methods of extracting information from poststenotic velocity disturbances and illustrates some of the advantages of applying modern signal analysis methods to fluid dynamic measurements.

REFERENCES AKAIKE, H. A. (1969). Fitting autoregressive models for prediction. Ann. Inst. Stat. Math. 21, 243. BEAVERS, C. S. ,g" WILSON, T. A. (1970). Vortex growth in jets. J. Fluid Mech. 44, 97. BERGLAND, G. D. (1969). A guided tour of the fast Fourier transform. IEEE Spectrum July 6, 41. BLACKMAN, R. B. & TUKEY, J. W. (1959). The measurement of power spectra from the point of view of communications engineering. New York: Dover. Box, G. E. & JENKINS, G. H. (1976). Time series analysis: forecasting and control. San Francisco: Holden Day. BURG, J. P. (1967). Maximum entropy spectral analysis. 37th Ann. Int. Meeting Soc. Explor Geophysics, Oldahoma City. Burg, J. P. (1972). The relation between maximum entropy spectra and maximum likelihood spectra. Geophysics 37, 375. CASSANOVA, R. A. & GIDDENS, D. P. (1978). Disorder distal to modeled stenoses in steady and pulsatile flow. J. Biomech. 11, 441. EDWARD, J. A. & FITELSON, M. M. (1973). Notes on maximum entropy processing. IEEE Trans. Information Theory IT-19, 232. GIDDENS, D. P. & KITNEY, R. L (1983). Autoregressive spectral estimation of poststenotic blood flow disturbances. J. Biomech. Eng. 105, 401-404. GIDDENS, D. P. MABON, R. F. & CASSANOVA, R. A. (1976). Measurements of disordered flows distal to subtotal vascular stenoses in the thoracic aortas of dogs. Circ. Res. 39, 112. GOSLING, R. G.) DUNBAR, G., KING, D. H., NEWMAN, D. L., SIDE, C. D., WOODCOCK) J. P., FITZGERALD, D. E, KEATES, J. S. & MACMILLAN, D. (1971). Angiology 22, 52. HAYKIN, S. S. (ed.) (1979). Nonlinear method of spectral analysis. New York: Springer-Verlag. KAY, S. M. (1979). The effects of noise on the autoregressive spectral estimator. IEEE Trans ASSP 5, 478. KAY, S. M. cg~MARPLE, S. L. Spectrum analysis--a modem perspective. Proc IEEE 69, 1380. KHALIFA, A. M. A. & GIDDENS, D. P. (1978). Analysis of disorder in pulsatile flows with application to poststenotic blood velocity measurement in dogs. J. Biomech. 11, 129. KHAL1FA, A. M. A. 8/. GIDDENS, D. P. (1981). Characterization and evolution of poststenotic flow disturbances. J. Biomech. 14, 279. KITNEY, R. I. & GIDDENS, D. P. (1982). Extraction and characterization of underlying velocity waveforms in poststenotic flow. IEEE Proc. Part A 129, 651. KITNEY, R. I. & GIDDENS, D. P. (1983). Analysis of blood velocity waveforms by phase shift averaging and autoregressive' spectral estimation. J. Biomech. Eng. 105, 398. KNOX, R. A., GREENE, F. M.) BEACH, K., PHILLIPS, D. J., CHIKOS, P. M. & STRANDNESS, D. E. JR (1982). Computer based classification of carotid arterial disease: prospective assessment. Stroke 13, 589. LINKENS, D. A. (1979). Empirical rules for the selection of parameters for autoregressive spectral analysis of biomedical rhythms. Signal Processing 1, 243. PAP.ZEN, E. (1968). Statistical spectral analysis (single channel case). In: Dept of Statistics--Stanford Univ Technical Rep 11. OPPENHEIM,A. V. & SCHAFER,R. W. (1975). Digital signal processing. New Jersey: Prentice Hall.

442

R. I. K I T N E Y E T AL.

ULRYCH, T. F. (1972). Maximum entropy power spectrum of truncated sinusoids. J. Geophysical Res. 77, 8, 1396. UUt¥CH, T. F. & BISHOP, T. N. (1975). Maximum entropy spectral analysis and autoregressive decomposition, Rev. Geophys. Space Phys. 13, 183. ULRYCH, T. J. & CLAYTON, R. W. (1976). Time series modelling and maximum entropy. Phys. Earth Planetary Interiors 12, 188.