Technical Field

This invention pertains to the measurement of blood flow by the analysis of Doppler shifted ultrasonic signals. In particular, the invention relates to analyzing ultrasonic signals utilizing autoregressive techniques in conjunction with synchronizing electrocardiogram (EKG) signals by averaging reflection coefficients resulting from the autoregressive analysis over a number of cardiac cycles. Background of the Invention

Due to the large number of deaths and physical disabilities arising from diseases of the heart and major blood vessels in both adults and children, there exists a well-recognized need for reliable, accurate, and inexpensive blood flow measurement techniques. There are two general types of methods fo-r measuring blood flow in humans - invasive and non-invasive techniques. - Since the invasive methods involve radioactivity and injections, only one study is performed per patient, and this method is not suitable for utilization as a diagnostic screening technique. There are two types of non-invasive blood flow measurements: Doppler flow meter and plethysmography. The plethysmography method involves the compression of the neck veins and is of little practical use for diagnostic screening.

The Doppler method is the most commonly used non-invasive method for performing blood flow analysis and can be repeated as necessary. In general, this method involves the transmission of an ultrasound signal through the skin to a blood vessel and the detection of

_ .

the Doppler shift in the reflected ultrasound signal resulting from the movement of the red blood cells. The Doppler shifted ultrasound signal is then utilized to determine the velocity of blood flow. The latter determination is complicated due to the presence of noise in the Doppler shifted signal, and because there are sets of blood cells moving at different velocities. Each of these sets of blood cells gives rise to a different frequency shift resulting in a Doppler shifted signal that is a complex wave. One example of noise that is present in the Doppler shifted signal is that due to the movement of blood vessel walls. The wall of a blood vessel moves out during the systole portion and returns during the systole portion of the cardiac cycle. These movements result in both low- and high-frequency noise components in the Doppler shifted signal.

Yet another problem that arises in attempting to determine the amount of blood flow from the Doppler shifted signals is that the heart rate varies not only from person to person, but within a given person over a relatively small number of cardiac cycles. This variation makes the interpretation of the resulting Doppler shifted signals difficult since the difference in cardiac cycles must be included if the blood flow information is to be extracted from the Doppler shifted signal since the resulting display from these signals can have a very different appearance.

Because of the complexity of the Doppler shifted signal, the presence of noise in the Doppler shifted signal, and the variation in cardiac cycles, it is necessary to perform analog or digital analysis of these signals in order to present the information in a manner which is useful to medical personnel attempting to diagnose blood flow rate. In small medical multiphase screening laboratories, the requirement of either having to immediately process and average data over a large number

of cycles and plot the resulting power spectra on a hard copy device or to store all of the resulting power spectras requires computer power and peripheral devices not always available in those type of laboratories. The problem is that the technique utilizing autoregressive analysis requires the calculation of the power spectra for each channel of every cycle and the large amount of data on an individual patient basis requires a large amount of storage. Summary of the Invention

In an illustrative method and structural embodiment, a departure in the art is achieved by performing an autoregressive analysis of reflected Doppler shifted ultrasonic signals resulting from the blood cell movement to determine the amount of flow by averaging the reflection coefficients and residual energy levels resulting from the autoregressive analysis. Each cardiac cycle of a patient is divided into a predefined number of time segments or channels by first determining the start and the end of each cardiac cycle from the patient's electrocardiogram signals. An autoregressive analysis is then performed on each individual channel to determine 'the reflection coefficients and the residual energy level for each channel. The reflection coefficients resulting from the autoregressive analysis closely model the blood flow within the patient's internal blood vessel and can be performed in an economical manner. The residual energy represents energy not accounted for by that modeling. The reflection coefficients and residual energy level for individual channels are then averaged over all of the cycles resulting in a small amount of data that represents the blood flow within the patient's blood vessel. These averaged reflection coefficients and residual energy levels can then be stored in economical mass storage devices for later analysis and display. For interpretation by medical personnel, the reflection

coefficients and residual energy level data are utilized to calculate averaged power spectra for individual averaged channels. Advantageously, these averaged power spectra are then displayed on a channel basis using different colors to represent different power levels thus greatly enhancing the ease of interpretation by medical personnel.

Advantageously, an autoregressive analysis system for obtaining the average power spectrum of a plurality of -signal groups each comprising a plurality of information signals functions by a first set of steps determining the start and end of each of set signal groups, and a second set of steps is responsive to the information signals of each group for calculating reflection coefficients. A third set of steps is responsive to the information signals and the reflection coefficients for a given group to generate a residual energy level for each group. After all of the reflection coefficients and residual energy levels for all the groups have been calculated, a fourth set of steps computes the average of the reflection coefficients and the average of the residual energy levels over all the groups. A fifth set of steps is responsive to the averaged reflection coefficients and averaged residual energy level for generating an average power spectrum.

Advantageously, the autoregressive analysis system utilizes a Doppler shifted flow meter to transmit a beam of ultrasonic frequency acoustical waves into an internal patient vessel and for receiving the reflected acoustical echoed signal back from moving blood cells within the internal patient vessel and determining the difference between the transmitted and received signals. An analog-to-digital (ADC) converter is used to convert the sampled difference signals into digitized signals which are the information signals of the signal groups.

In addition, an electrocardiogram instrument is used to obtain the patient's electrocardiogram signals, and a second ADC converter digitizes the electrocardiogram signals. The first set of steps is responsive to the digitized electrocardiogram signals for determining the start and end of the signal groups which correspond to the start and end of cardiac cycles. Advantageously, the second set of steps which calculates the reflection coefficients comprises three subsets of program instructions. The first subset of steps computes a forward error between each of the information signals of each signal group and subsequent ones of the information signals of the same group, and a second subset of steps computes backward error signals for each of the information signals of each group representing the difference between each information signal and the preceding information signals of the same group. Finally, the third subset of steps computes the reflection coefficients in response to the ^" backward and forward error signals for each of the groups.

A method for determining the amount of blood flow by utilization of a system comprising a Doppler flow meter, electrocardiogram instrument, analog-to- digital converters, and a digital computer performs the following steps: determines the start and end of a patient's cardiac cycles in response to electrocardiogram signals by an analysis of digitized signals from the electrocardiogram instrument, digitizing Doppler shifted signals from the flow meter by utilization of an analog-to-digital converter, dividing the resulting digitized Doppler shifted signals into a predefined number of time channels, calculating reflection coefficients for each of said channels for said Doppler shifted signals, computing residual energy level in response to the calculated reflection coefficients and the Doppler shifted signals for each of said channels, averaging the reflection coefficients and

the residual energy levels for particular channels over all of the cardiac cycles, and calculating the power spectrum from the averaged reflection coefficients and averaged residual energy level for each channel. Advantageously, the method further comprises the steps of normalizing the average power spectrum for each of the combined channels and displaying the normalized power spectrum for each channel utilizing different colors to represent f e different power levels thus enhancing the display for the ease of medical personnel performing an analysis of the amount of blood flow.

Brief Description of the Drawing

FIG. 1 illustrates, in block diagram form, a blood flow analysis system in accordance with this invention;

FIG. 2 illustrates in graphic form, the results of analysis performed by the blood flow system illustrated in FIG. 1 by averaging the .reflection coefficients and residual energy that model the flow of blood within an artery of a piglet;

FIG. 3 illustrates, in flowchart form, the steps executed by computer 106 of FIG. 1 in performing the blood flow analysis; FIG. 4 is a detailed flowchart of the steps executed by blocks 302 through 304 of FIG. 3;

FIGS. 5 and 6 illustrate in greater detail, the steps of block 409 of FIG. 4;

FIG. 7 is a detailed flowchart of blocks 307 through 310;

FIG. 8 is a detailed flowchart of block 708 of FIG. 7; and

FIG. 9 is a detailed flowchart of block 709 of FIG. 7.

Detailed Description

A system for analyzing and displaying the flow of blood within a blood vessel is illustrated in FIG. 1. Doppler flow meter 102 and ultrasonic transducer 101 generate an analog Doppler shifted signal that corresponds to the speed of the red blood cells (RBC) that are moving in a blood vessel. The frequency of the signal indicates the speed at which the blood cells are moving, and the power -cr- served at a particular frequency indicates the number of blood cells moving at that frequency. Because of such factors as the angle of the probe and other complex physical factors, the energy does not literally correspond to the number of red cells. However, the energy is roughly proportional to the number of cells moving at any given speed.

The Doppler shifted signal is converted to digital form by analog-to-digital converter 103. Initially, computer 106 is responsive to the digitized Doppler shifted signal to store the signals along with the digitized EKG signal from electrocardiogram unit 104. After the information for a number of cardiac cycles has been stored, computer 106 is responsive to the digitized Doppler shifted signals and the digitized EKG signals to calculate the spectral power present at the various frequencies, to average the spectral ".powers for a number of cycles together, and to display this average spectral power over the time required for one cycle.

FIG. 2 illustrates the information that is displayed by computer 106 in analyzing the blood flow. Illustrated in FIG. 2 are the results of two tests performed on a piglet weighing approximately 4.2 kilograms and whose blood flow was artificially controlled. The blood flow rate was also monitored using electromagnetic flow meter that required surgical procedures and the results of the electromagnetic flow meter (EMF) are shown in the upper portion of each cycle

of FIG. 2. The difference between a flow rate of 50 milliliters per minute and a flow rate of 20 milliliters per minute, as determined by the EMF meter, is quite apparent in FIG. 2 and readily discernible by medical personnel. Each test is divided into a number of time segments referred to as channels; and illustratively, there by be 50 channels per cardiac cycle.

The display illustrates the average of the power on a per-channel basis for a plurality of cycles by first determining the reflection coefficients and the residual energy by utilizing autoregressive analysis techniques for each cycle, averaging the reflection coefficients and residual energy for individual channels for all of the cycles, and then determining the power level for each individual channel from the averaged reflection coefficients and averaged residual energy. The display illustrated in FIG. 2 is produced by accurately determining the start of each cardiac cycle and the end of the cycle by analyzing the EKG signals from the patient. Once the start and end of the cardiac cycle is determined, then the digitized Doppler shifted signals in this time period are illustratively divided up into 50 channels. An autoregressive analysis is performed on all of the Doppler shifted signals in a given channel to determine the reflection coefficients and residual energy of that channel. Some implementations of autoregressive analysis are the maximum entropy method of Berg or the Durbin Levinson Recursions on Autocorrelation Coefficients method. After all of the channels for one cycle have been processed, then the next cycle is processed. When all of the reflection coefficients and residual energies for all of the channels of all of the cardiac cycles have been determined, the average of the reflection coefficients and the residual energies of individual channels of all of the cycles is taken. Finally, the

power spectrum for each individual averaged channel is determined from the averaged reflection coefficients and residual energies for that channel and the resulting power spectra of the average channels is then displayed with different colors representing different power levels.

FIG. 3 illustrates in greater detail the steps necessary to produce a display such as illustrated in FIG. 2. For each cardiac cycle, block 301 determines the start and end of the cycle and the number of samples per channel. Reflection coefficients and residual energy for each channel are then determined by blocks 302, 303, and 304. Block 302 performs a autoregressive analysis of the samples for one channel to determine the reflection coefficients and the residual energy. The reflection coefficients define a model constructed by the autoregressive analysis of red blood cell movement. The residual energy is the amount of energy unaccounted for by the model. The reflection coefficients and residual energy are later averaged by blocks 305 and 306 and these averaged coefficients and energies are used to determine the spectrum for each averaged channel by block 307.

After the reflection coefficients and the residual energies have been calculated for all channels of all of the cycles, the average of the reflection coefficients and the residual energy for each particular channel of all of the channels is then calculated by blocks 305 and 306. The result of these averaging calculations is reflection coefficients and energy for each channel which represent the average for that channel over all of the cycles. Once the average of the reflection coefficients and the residual energies for the channels has been computed, the power spectrum of each channel is calculated for a multitude of frequencies utilizing the averaged reflection coefficients and residual energy for this particular

channel as calculated in block 305. For each cycle, block 307 and 308 determine the average spectrum for each channel. Once the average spectra for all of the cycles has heen calculated, the latter is displayed by block 311 on a color printer. Before display, the average spectra is normalized to take out variations due to the differences of the gain in the Doppler flow meter 102.

Autoregressive spectrum analysis of block 302 assumes an analysis model given by:

(1) E(z) = X(z) F(z) where: E(z) is the z-transform of the input signal,

X(z) is the z-transform of the output signal, and F(z) is an all pole filter.

In the sample domain, equation 1 can be written as follows :

(2)

ORDER e(n) = fc(m) x(n - m) m=0

which can be rewritten as

(3)

ORDER e(n) = x(n) - ^* " fc (m) x(n - m) m=l

where x(n) represents the present time sample, fc represents filter coefficients, and ORDER represents the number of filter elements.

As described in J. D. Markel and A. H. Gray, "Linear Prediction of Speech", Springer-Verlag, Berlin Herdelberg New York, 1980, on page 10, the input signal e(n) can be interpreted as the prediction error between the actual data sample x(n) and the linear combination

M of the previous n samples given by _ fc(m) x(n - m) . m=l Since the actual data sample minus the predicted signal can be interpreted as an error, the common analysis

procedure is to minimize the sum of the squares of this error as a method for determining the filter coefficients.

Many autoregressive techniques have been developed to minimize this error term; and the one utilized here is the Berg Maximum Entropy method which is described in the paper by L. Marple, "A New Autoregressive Spectrum Analysis Algorithm", IEEE Trans, on Acoustics, Speech, and Signal Processing, Vol. ASSP- 28, No. 4, August, 1980, pp. 441-454. The Berg method solves for the filter coefficients, fc(m) , by using a forward and backward prediction errors. For any given point, the er-ror is calculated by considering the points in time preceding the sample point under calculations and points prior to the particular point under consideration. The forward error is defined by:

(4)

ORDER f _prr [k) = S ^f c(i) * x(k - ORDER - i) ^{e r} i=0

while the backward error is defined by: (5)

ORDER b _{{k) =} _ fc(i) * x(k + i) i=0 where the values of index k range from 1 to (NPTS-ORDER) and fc(0) is defined as 1. NPTS is the number of sample points in the channel being analyzed. Using these two expressions, the problem is to determine the set of filter coefficients that minimize the sum of both the forward and backward errors summed over all sample points subject to a constraint. Thus, the error to be minimized is given by the following equation: (6)

NPTS-ORDER ^error = I ^f err- ^k > * f _e rr< ^k > ^{+ b} err< ^k > * ^b err ^(k » -

The constraint on the values of fc is that they satisfy the relation

(7)

fc(k) = FC(k) + re (ORDER) * FC(ORDER-k) where: fc(ORDER) = rc(ORDER) . where FC represents the coefficients determined when

ORDER-1 terms were used in the prediction equation.

This relation is called the Levinson recursion, and the term, rc(ORDER) , is frequently referred to as the LPC reflection coefficient.

The filter coefficients fc(k) can be solved for in equation 6 by substituting equations 4 and 5 into equation 6 and taking the partial derivative of this equation with respect to fc(ORDER) . The resulting equation is a recursive formula given in terms of the reflection coefficient as follows:

(8) NPTS-ORDER - ² . ^b _err (k)f _err (k+l) re (ORDER) =-

NPTS-ORDER ϊ (k+1) k=l err (k) ^* err

FIGS. 5 and 6 illustrate steps fo ^' r evaluating equation 8.

Given the filter coefficients, the power spectrum for any frequency, f, is given by

(9) P(f) = PE

ABS ₁₊ ° ^R 5 ^ER fc ₍ k ₎ *e(-2"* ⁱ * */fsample) k=l where fsample is the frequency at which the original signal was sampled, P(f) is the power at the particular frequency f, and PE is a measure of the residual energy left over in the autoregressive analysis. ABS represents the magnitude function of a complex number.

The steps illustrated in later material by flowcharts and program listings require equation 9 to be transformed to equation 10 for the following reasons. The power is computed in db to avoid large numbers. Thus, log(P(f)) is the quantity calculated. The exponential factor is broken down into real and imaginary parts. For these reasons, equation 9 is

transformed into the following:

(10)

10*log(P(f) )=10*log(PE)-20*log(ABS ORDER ( )

1+ £ fc(k)*|cos(k*f*arg)+i*sin(k*f*arg) I k=l ( ) where arg = 2 * π/fsample.

The equations 4 through 10 give a mathematical description of the process required to implement the flowchart of FIG. 3.

Now consider in greater detail, the various steps that are required to perform the previous equations. The process of performing the autoregressive analysis, as defined in block 302 of FIG. 3, is shown in greater detail in the flowchart illustrated in FIG. 4 as steps ulpc.c. As illustrated in FIG. 4, the steps first open the necessary files, block 401, in order to gain access to the digitized Doppler signal and the points which define the start and end of the cardiac cycles. Block 402 first reads the beginning of the first cycle and blocks 403 and 404 checks if the entire file has been read. If the entire file has been read, then the steps are finished. However, if the entire file has not been read, block 405 utilizes the starting point of the cardiac cycle to seek to the proper place in the Doppler signal file. Block 406 then computes the number of sample points for each channel by dividing the total number of samples by the number of channels. Illustratively, the number of channels is 50. Blocks 407 through 414 then compute the reflection coefficients and residual energy for a particular channel. For each channel, the data samples of that channel are utilized to calculate the reflection coefficients and the residual energy in block 409. The latter block implements the Berg method by utilizing fbeig steps which are illustrated in FIGS. 5 and 6. The reflection coefficients are then stored in an output file by block 410. The reason for storing reflection coefficients rather than the filter coefficients

themselves is that the reflection coefficients are bounded between 1 and -1 whereas the filter coefficients can be quite large or small. The residual energy is converted to db and stored by blocks 411 and 412. Blocks 413 and 414 determine whether or not all the points of the channel have been utilized. If they have, then the steps return to block 403 to process another cycle.

The fberg steps is illustrated in greater detail in FIGS. 5 and 6. These steps calculate equation 8 and verifies that the conditions of equation 7 have been met. The fberg steps calculate equation 8 in a recursive manner by first calculating the first reflection filter coefficient for order 1 and then uses this information to calculate a new set of filter coefficients for order 2. Because of this recursive nature, fberg steps are repetitively recalculating equation 8. Blocks 501 through 506, illustrated on FIG. 5, perform the initialization of the various variables used by the subroutine. The initial residual energy, etot, is set equal to the sum of the squares of all the data points of the channel and the backward and forward errors (berr and ferr, respectively) are set equal to corresponding data samples, where the data sample is the digitized Doppler signal. The filter coefficients are initially set equal to 1 by block 506 and the denominator of equation 8 is set equal to the initial residual energy as determined by block 503. FIGS. 5 and 6 illustrate the evaluation of equation 8. Block 507 determines whether or not the equation has been sufficiently evaluated for the order filter being calculated, and if it has, a return is executed to the program illustrated on FIG. 4 via return block 508. If the filter has not yet been calculated to a sufficient order, then block 509 is executed. The numerator and denominator (num and den, respectively) of

equation 8 for this particular order are evaluated by blocks 510 through 513. After the numerator and denominator have been determined, then the reflection coefficient, rc[m] , and residual energy for this particular order are evaluated in block 514. Once the reflection coefficient for this order has been determined, then the stepup function is implemented by blocks 515 through 519 to update the previously determined filter coefficients for (ORDER -1) as defined by the evinson recursion formula and given in equation 7. The highest order filter coefficient is always equal to the reflection coefficient and is set equal to the reflection coefficient by block 520. The past filter coefficients are then updated by blocks 521 through 523 in order for block 517 to evaluate the next set of filter coefficients. The past filter coefficients are designated as FC(n) in equation 7, and as pfc[n] on FIGS. 5 and 6. The forward and backward errors are next updated by blocks 524 through 526 in order to evaluate blocks 511 and 513 in the next iteration. After the forward and backward errors have been updated, control is passed from block 526 to decision block 507 which determines whether or not all the orders have been evaluated. Referring back to FIG. 3, after the reflection coefficients and residual energy levels have been calculated in blocks 302 through 304, the average for each individual channel over all cycles is performed as illustrated in blocks 305 and 306. Again, referring back to FIG. 3, once the averaging has been performed by blocks 305 through 306, the power spectrum for each channel is calculated for a plurality of frequencies as defined in equation 10. Illustratively, for each channel, equation 10 is evaluated for 150 different frequencies by blocks 307 and 308. After the average power spectra has been calculated for the cycles, the resulting power spectra is normalized by block 310 and

then displayed by block 311.

Blocks 305 through 310 are implemented by the cfdspect.c steps which are illustrated as a flowchart in FIG. 7. Blocks 701 through 704 perform the initialization required for the execution of the cfdspect.c program. The sine and cosine values of equation 10 are initially evaluated and stored so as to save processing time. The av array which is used to hold the averages for the spectra for the different channels is initialized to 0. Next, the arrays which hold the values for the reflection coefficients and the residual energy levels for all the channels of all the cycles are initialized by block 715 reading this information from a file created by blocks 410 and 412 of FIG. 4. Next, block 716 sums the reflection coefficients and residual energy level for each particular channel for all cycles with reflection coefficients for the same filter element being summed together for all cycles. After the reflection coefficients and residual energy levels have been summed for all cycles, the power spectrum for each combined channel is calculated by blocks 706 through 710 in the following manner. For each combined channel, the average of the reflection coefficients and the residual energy level is calculated by block 707. Block 708 then converts the reflection coefficients for the channel into filter coefficients by executing the stepup subroutine that implements equation 7. These steps are illustrated in FIG. 8. Once the filter coefficients have been calculated, the power spectrum is calculated by block 709 which uses the cspect steps of FIG. 10. Blocks 706 and 706 control the steps so that the previous calculations are performed for all of the combined channels. Once the spectra for all of the combined channels has been calculated, then blocks 711, 728, and 729 are executed implementing block 310 of FIG. 3

which normalizes and stores the resulting spectra of all of the cycles. The normalized spectra are displayed on a color printer.

It is to be understood that the above- described embodiment is merely illustrative of the principles of this invention; other arrangements may be devised by those skilled in the art without departing from the spirit and the scope of the invention.