Application of DSP Concept for Ultrasound Doppler Image Processing System

Abstract

Blood-flow measurements using Doppler ultrasound system are popular in ultrasonic diagnoses. But the blood-flow measurement inside the heart is difficult. There are many reasons behind it. The deep range and fast blood-flow are difficult to measure because of limitation of acoustic velocity. Moreover, strong heart valve signals mix into the blood-flow signal. Against such difficulties, the statistics mathematical model was applied to analyze many clinical data sets. The system identification method based on the mathematical model could realize a new blood-flow measurement system that has ultrasound Doppler information as input and electrocardiogram as output.

1. Introduction

It has been more than 70 years since Doppler ultrasound technology was born [1]. In the meantime, in the field of medical blood-flow measurement, various diagnostic methods and diagnostic indices were proposed and had been standardized. By re-focusing on the empirical data from these diagnoses, the relationship between biosignals that many medical doctors had accumulated was revealed. In this chapter, new applications of statistical diagnosis methods and imaging technology are proposed.

Medical Doppler ultrasound systems are commonly used for various diagnostic applications, including examination of cardiac and abdomen. Figure 1 is the example of diagnostic image of a left ventricle inflow. The upper part of Fig. 1 (B-mode image) shows the tomogram of echo, and the lower part of Fig. 1 (D-mode image) shows the spectrum Doppler image. The D-mode image shows the blood-flow velocity at the mitral valve on the B-mode image. In the D-mode image, the horizontal axis is time and the vertical axis is the blood-flow velocity which corresponds to Doppler-shift frequency. The waveform displayed in the lower part of D-mode image is an electrocardiogram (ECG). The amplitude of echo reflected from tissue constructs B-mode image, and the Doppler-shift signal from blood-flow constructs D-mode image.

Figure 1.

Diagnostic image of Doppler ultrasound system.

2. A new approach to the Doppler diagnostic: ARX model and biosignals

2.1. Accumulated medical database

Currently, time-sharing blood-flow measurement in Doppler ultrasound system appeared. But the time-sharing systems have many problems caused by acoustic velocity range limitation. To address these problems, the mathematical models based on system identification methods were proposed in this chapter. One of the system identification models has ECG as input and has short time Fourier transform (STFT) image parameters as outputs. Based on this model, a new gap-filling system introduced in Section 3 was developed. It can fill the 100 ms gap.

Doppler ultrasound diagnoses of the left ventricle inflow and outflow are very helpful. The diagnostic techniques using the peak velocity waveform (the envelope trace of the Doppler spectra) are the standards of blood-flow measurement. Synchronizing with the systolic phase and diastolic phase that ECG shows, the mimetic diagrams of the outflow from an aortic valve and the inflow from a mitral valve are shown in Fig. 2. The features of these waveforms are measured and evaluated, and they are used as the standards of cardiac disease diagnoses [2]. New diagnostic technology that applies causal relationship between biosignals (here, they are an ECG and a Doppler trace waveform) introduced in Section 4 was developed. Furthermore, many medical doctors make standards of causal relationship between these biosignals over the time of 30 years or more. They are suitable to be applied to the statistical models.

Figure 2.

Heartbeat indices and their measurement guidance.

2.2. Mathematical models used in system identification

In order to express the causal relationship between biosignals, the mathematical model that consists of an input, an output, and a black-box is suitable. The black-box model is shown in Fig. 3. Since the input x and the output y assume multi-inputs and multi-outputs (MIMO), they are expressed as vectors. Moreover, in order to take a time response into consideration, two linear partial differentiation equations (1) and (2) express the model using the state variable u.

Figure 3.

Black-box model and linear dynamic system.

$$\dot{\overrightarrow{u}}(t)=\overrightarrow{A}\cdot \overrightarrow{u}(t)+\overrightarrow{B}\cdot \overrightarrow{x}(t)$$E1

$$\overrightarrow{y}(t)=\overrightarrow{C}\cdot \overrightarrow{u}(t)+\overrightarrow{D}\cdot \overrightarrow{x}(t)$$E2

ECG, the Doppler waveform and the spectrum Doppler image were used for the causal relationship analyses. In Section 3, the mathematical model based on ECG and Doppler imaging parameters was used. In Section 4, the mathematical model based on ECG and Doppler trace waveform was used. Figure 4 shows the expression of ARX model frequently used in these investigations.

Figure 4.

Numerical formula of ARX model.

2.3. System identification based on ECG and Doppler waveform

When ECG and the Doppler waveform are applied to ARX model, the system is expressed as the model of Fig. 5. Two coefficient sequences a_i and b_j which determine the system response are calculated using a statistical method. This is the system identification using ARX model.

Figure 5.

ARX model based on ECG and the Doppler waveform.

The physical meaning of ARX model can be explained with IIR type digital filter. Figure 6 is a digital filter with two inputs: noise(n) and ECG. Since noise(n) is random, there is no time causal relationship, but ECG has FIR ingredients. V_p (Doppler waveform) has IIR ingredients of its feedbacks.

By using system identification, the causal relationship between ECG and V_p was summarized in two coefficient sequences. These coefficient sequences are equivalent to the feedback coefficient sequence and feed-forward coefficient sequence of the digital filter. The waveform of ECG resulting from the pulsation of a heart can be explained by how it is related to the flow velocity expressed by the Doppler waveform. It is also the same expression as IIR digital filter.

Thus the response of black-box model can be presumed by system identification using statistical data. Noise rejection technology of the Doppler waveform and gap-filling technology of Doppler image based on system identification are introduced. Moreover, possibilities such as technology that complements the lack part of ECG and automatic diagnostic technology (computer aided diagnosis [CAD]) using statistical data will be introduced to another opportunity.

Figure 6.

ARX model denoted by digital filter expression.

3. Application of system identification to a new gap-filling algorithm

3.1. Limitation of ultrasound scans by acoustic velocity

Doppler velocity range limitation caused by transmit pulse repetition frequency (PRF) is a serious matter with blood-flow measurement. There exists a trade-off between depth range and velocity range.

In order to display the B-mode image and the D-mode image simultaneously in real time, a time-sharing scanning is needed. Normally, PRF of approximately 4 kHz is employed, taking the propagation time of acoustic wave and the attenuation in the living body into account. Use of higher PRF has many advantages. For example, the scanning line density is increased, and as a result B-mode images with higher azimuth resolution can be obtained. In addition, the velocity range of D-mode is expanded. On the other hand, the higher PRF reduces the imaging depth range. Therefore, information concerning deeper regions cannot be obtained. So PRF control is complicated, especially with time-sharing of B-mode scanning and D-mode scanning.

In current Doppler ultrasound systems, it is difficult to optimize both the D-mode image quality and the B-mode image quality simultaneously. A new Doppler gap-filling algorithm based on ARX model was developed, which had ECG as an external input, for detecting the high-speed blood-flow in heart, carotid arteries, etc., and for generating high-quality D-mode images. The conventional gap-filling algorithm of D-mode image suffers from various problems such as presence of noise or artifacts and poor reproducibility in the rapid velocity change.

Figure 7.

Examples of the interleave scanning and the segment scanning.

Examples of time-sharing transmission/reception of the interleave scanning and the segment scanning are shown in Fig. 7. The B-mode images (100 beams, 6 kHz PRF, approximately 7 cm depth) with a frame rate of 30 Hz are obtained in both Fig. 7(b) and (c). However, the velocity ranges obtained simultaneously in D-mode differ. The velocity range of D-mode is set to 6 kHz in the interleave scanning shown in Fig. 7(b). But the velocity range of D-mode is set to 12 kHz in the segment scanning shown in Fig. 7(c). Since both B-mode and D-mode become discontinuous in the case of segment scanning, the gap-filling algorithm is needed in D-mode signal processing and interpolation processing is needed in D-mode image processing. Moreover, when the gap-filling algorithm is applied, both quality of the image and the audio are degraded. On the other hand, since the PRFs can be set to B-mode and D-mode independently, the Doppler velocity range can be expanded.

3.2. Conventional gap-filling algorithm

Simultaneous real-time display of B-mode and D-mode by the segment scanning has been used for many years. The gap-filling algorithm fills in the gaps of IQ signal (shown in Fig. 8). The gap is filled with the predicted waveform that is generated based on an autoregressive (AR) model. Recently, many improved gap-filling algorithms have been reported [3]. For example, in one method, the gaps are filled in from both time directions; in another method, narrow-band noise is used as a source of signal; and in another method, an autoregressive moving average (ARMA) model is used.

Figure.8 shows a Doppler ultrasound system with the conventional gap-filling algorithm. The received beam is generated in digital beam former (DBF). The output of DBF is sent to the envelope detection processing in which echo intensity is detected. Then the echo intensity signal is sent to the B-mode image processing, and then displayed as the B-mode image. In the spectrum Doppler processing section, the Doppler-shift signal is detected by quadrature detection. Since the detected signal contains low-velocity and high-power components called clutter from tissues (vessel walls etc.), the wall filter rejects the clutters except for blood-flow components. The gap-filling algorithm interpolates the gaps of the D-mode image and the Doppler audio caused by the segment scanning.

Figure 8.

Doppler ultrasound system with conventional gap-filling algorithm.

Figure 9 shows the details of gap-filling algorithm. It shows system identification and linear prediction based on AR model. Figure 10 shows the timing chart of its signal-processing shown in Fig. 9. To estimate the output for Gap(B), band limiting is applied to a white noise source. Prediction is performed in both forward and reverse directions, and blending is performed with W₁(t) and W₂(t) in order to improve continuity. W₁(t) and W₂(t) are weighing functions used for blending actual waveform and predicted waveform. Predicted waveforms of two directions fill in the Gap(B). Because the rapid audio response is important, generally only forward prediction is used. A part of the time sequential data before the gap is used to calculate forward prediction coefficient series Af(p) and bandwidth (BW) of the residual noise. Using Af(p) and BW obtained from the data immediately before Gap(B), the IQ signals are estimated and outputted.

Figure 9.

System identification and prediction based on AR model.

Figure 10.

Timing chart of gap-filling algorithm.

3.3. Problems of conventional system

There are two problems specific to Doppler ultrasound systems employing the conventional gap-filling algorithm. The first problem is that artifacts are newly generated in the predicted output although the low-frequency components already have been removed in the wall filter. For these artifacts, not only D-mode image quality but also audio quality is degraded. The second problem is that discontinuity of D-mode images becomes greater when there are sudden changes in the spectrum (rapid changes in blood-flow velocity). In many cases, the predicted D-mode image has horizontal lines and spikelike noises that are observed near the gaps. Figure 11 shows the D-mode images of a portal vein with moderate changes in velocity. Figure 11(a) is the D-mode image obtained by interleave scanning and Fig. 11(b) is the D-mode image obtained by segment scanning with the conventional gap-filling algorithm. Fig. 11(a) is smooth and free of artifacts, while periodic spikelike noise and low-frequency components are observed in Fig. 11(b).

Figure 11.

Artifacts in D-mode image caused by the conventional gap-filling algorithm.

The conventional gap-filling algorithm based on AR model (or ARMA model) uses the noise and the predicted output itself as feedback inputs. So it tends to generate the waveforms consisting of multiple changeless frequency components. When noise is x(n), output is y(n), and coefficient series obtained in system identification is a_k, the predicted output is expressed by equation (3). In the gap-filling algorithm based on AR model, the noise with Gaussian distribution or the noise with narrow bandwidth is used for n(n). Assuming that σ² is distribution width of noise n(n), the estimated spectrum output P(f) is expressed by equation (4) [4]. T is the sampling time. Equation (4) suggests that AR model can only generate multiple spectral peaks in steady states. This is not suitable for estimating rapid changes in velocity. It has been reported that the time to be considered as steady states in D-mode is approximately 10 ms. Accordingly, conventional systems have been designed to limit the Gap(B) to 10 ms or less. It is also known that if the Gap(B) is longer, number of beams and continuity in B-mode image increase and image quality is improved.

$y\left(n\right)={\displaystyle \sum _{k=1}^P{a}_k\cdot y\left(n-k\right)}+n\left(n\right)$E3

$P\left(f\right)=\frac{{\sigma }^2\cdot T}{{\left|1-{\displaystyle \sum _{k=1}^P{a}_k\cdot \exp \left(-j\cdot 2\pi \cdot f\cdot T\left(k\right)\right)}\right|}^2}$E4

To improve image quality of B-mode, a longer gap than D-mode is required. But D-mode image quality is markedly degraded when the gap of B-mode becomes long. It is also important to track blood-flow changes due to pulsation for D-mode image quality. But the gap-filling algorithm based on AR model is insufficient. Diagnostic performance will be substantially improved if long gaps are not filled with changeless spectra but filled with changeful spectra.

3.4. A new gap-filling algorithm based on ARX model

A new algorithm that can reduce spectrum artifacts and stabilize rapid changes in velocity in order to overcome problems shown in Section 3.3 was developed. This algorithm uses not IQ signals but Doppler spectrum parameters as input, and is based on ARX model [5]. The outline of a new D-mode image processing is shown in Fig. 12. After quadrature detection IQ signals are generated. IQ signals are processed by the wall filter and STFT sequentially, and D-mode image is generated. The waveforms with 600 ms time lack (left time-domain IQ signals) show the output of the wall filter, which removes low-frequency clutter. The output of STFT shows the momentary spectra (right frequency-domain periodgram). STFT conducts frequency analysis and carries out the time shift image of spectra [6, 7].

Figure 12.

New D-mode image processing.

A new gap-filling algorithm in Fig. 12, which is based on ARX model, is shown in Figs. 13 and 14. Figure 13 is a block diagram of system identification for the new gap-filling algorithm. Figure 14 shows ECG and the D-mode image of left ventricular inflow. The spectrum shown in the lower side of Fig. 14 shows mean velocity (V_m) and distribution (σ) and the model spectrum near the time of 1.5 s in D-mode image. In system identification based on ARX model, time sequence data (coefficient series) during the cardiac cycle is calculated. First the power spectrum SP(f,t) is calculated, and then V_m, σ, and spectrum total power (TP) is calculated. The V_m, σ, and TP are the characteristic parameters of a single-peak spectrum model this time.

Figure 13.

System identification of a new algorithm.

Figure 14.

Single-peak spectrum parameter model.

3.5. System design and its confirmation

The formulas of spectrum parameters are shown in equations (5)-(7). For example, the power spectrum is expressed as P(f) and number of FFT points is set to 128. Here, f is frequency, and FFT sampling frequency f_s(=1/T) is 128. The value of f therefore ranges from 0 to 127*f_s/128. Left ventricular inflow shown in Fig. 14 is used as the input. Only positive side components are used in calculation, so value of k ranges from 0 to 95.

$V_m=\left(\frac{C}{2}\right)\cdot \frac{{\displaystyle \sum _{k=0}^{95}{f}_k\cdot P\left(f_k\right)}}{{\displaystyle \sum _{k=0}^{95}P\left(f_k\right)}}$E5

$\sigma =\sqrt{\frac{{\displaystyle \sum _{k=0}^{95}\left(f_k-\frac{2\cdot V_m}{c}\right)}}{96}}$E6

$TP={\displaystyle \sum _{k=0}^{95}P\left(f_k\right)}$E7

The system that has V_m(n), σ(n), and TP(n) waveforms as outputs and ECG waveform as an input was applied to ARX model. These parameters can be obtained by system identification shown in Fig. 13. At the same time these parameters are used for the new gap-filling algorithm. Figure 15 shows the block diagram of the spectra prediction processing based on ARX model. This time ECG was newly added as an external input, and only forward prediction was used in the mathematical model.

Figure 15.

New gap-filling algorithm based on ARX model.

The predicted output of ARX model $\overrightarrow{y}\left(n\right)$ is shown in equation (8). Here, a_k and b_k are ARX model coefficients. $u\left(n\right)$ is the normalized (both time and amplitude) ECG waveform.

$$\overrightarrow{y}(n)={\displaystyle \sum _{k=1}^P{a}_k}\cdot \overrightarrow{y}(n-k)\begin{array}{cc}& +\begin{array}{cc}{\displaystyle \sum _{k=1}^q{b}_k}\cdot u(n-k)\begin{array}{cc}& +\end{array}& n(n)\end{array}\end{array}$$E8

$$\overrightarrow{y}(n)$$ is a vector and is expressed by equation (9) as follows:

$$\overrightarrow{y}(n)=\left[\begin{array}{c}{V}_m\left(n\right)\\ \sigma \left(n\right)\\ TP\left(n\right)\end{array}\right]$$E9

The first term and third terms of equation (8) are similar as equation (3). However, the second term of equation (8) means a time-variant system and can generate changeful spectra.

Simulations were applied to confirm ARX model processing and its performance. The data used for simulation was left ventricular inflow, which exhibits rapid changes in velocity. The gap of segment scan was set to 100 ms, which is one order of magnitude larger than the conventional system. The simulation result is shown in Fig. 16. Time 0 to 1 s is a continuous D-mode image (without segment scanning), and time 1 to 2 s is a discontinuous D-mode image (segment scanning). Domains indicated by (a) are actual spectra, and domains indicated by (b) in the figure are spectra estimated by ARX model. Condition of the simulation is as follows: ARX prediction order is 9 to 12, the segment gap is 100 ms, and the blending time is 16.7 ms.

Figure 16.

Simulation of predicted spectra based on ARX model.

This result shows that spikelike noises and low-frequency artifacts are reduced compared with the conventional algorithm. Moreover, it is possible to obtain a changeful D-mode image under conditions of rapid changes in velocity and larger segment gaps.

4. Application of system identification to automatic cardiac valve-rejection algorithm

4.1. Problem of blood-flow measurements

The blood-flow diagnoses by Doppler ultrasound system have become popular recently. Peak velocity of blood-flow (STFT envelope waveform) is traced automatically in this system. But valve signals are mixed with the blood-flow signals in the heart. So automatic blood-flow measurements are not correctly recorded. To solve this problem, the mathematical model that has ECG as an input and has Doppler waveform as an output was applied. Using system identification method, a new valve-rejection algorithm was developed [8, 9].

Figure 17 is a Doppler ultrasound diagnostic image for left ventricular outflow. Doppler ultrasound system traces V_p (peak velocity of Doppler waveform) automatically. V_p is superimposed as a bright yellow line. In cardiac blood-flow measurements, cardiac wall noises with strong and low-velocity ingredient generate low-velocity artifacts. The valve noises with strong and high-velocity ingredient generate spikelike trace errors. Although cardiac wall noises seldom influence V_p waveform, valve noises have unneglective influence on automatic tracing. Valve noises usually mixed in left ventricular outflow in Fig. 17. Thus, users must compensate V_p waveform manually based on the Doppler spectrum image.

In Fig. 17, ECG (green line) is displayed with V_p waveform simultaneously. Also the R-triggers (sharp peaks of ECG) are at 1.48 s and 0.65 s. The systolic phase from the R-trigger of ECG is about 300 ms. Conventional manual trace (cyan line) is displayed between 1.5 s and 1.1 s. Medical doctors estimate left ventricular outflow based on manual traces.

Figure 17.

Example of left ventricular outflow.

4.2. Blood-flow measurement

Figure 18 shows a Doppler ultrasound system and automatic blood-flow measurement system. DBF generates echo beams from transmitted and received ultrasound signals. V_p waveform is the spectra envelope of D-mode image, and it is automatically traced and superimposed. B-mode image and D-mode image, and V_p and ECG are simultaneously displayed in the same screen. Various blood-flow measurements that combine V_p and ECG are known in clinical applications.

Figure 18.

Automatic blood-flow measurement system.

4.3. System identification of left ventricle outflow

Many clinical data sets (V_p, ECG, diagnostic indices and image) were acquired from numerous volunteers. Using these data sets, system identification shown in Fig. 19 was investigated. V_p has both valve and blood-flow information, so valve components were removed from V_p manually and the ideal blood-flow waveform (V_i) was generated. There were variations in these data sets, such as heartbeat cycles and flow velocity ranges. They were normalized for each data set. Amplitude of V_i was normalized by maximum flow velocity of the systolic phase. Amplitude of ECG was normalized by its R-trigger voltage. Heart-beat cycles were normalized and divided by 60. So sampling periods were fixed to 60 (approximately 60 Hz). A low-pass filter with a one-third cutoff (approximately 20 Hz) was applied after normalization, and unnecessary frequency components were rejected. A normalized V_i (NV_i) and a normalized ECG (NECG) were obtained after filtering process. The system identification block has NECG as an input and has NV_i as an output. The coefficient sequences of mathematical model were generated by system identification. A diagnostic image of left ventricular outflow is shown in Fig. 20. V_p and V_i (the aortic valve signal was rejected manually) are shown in Fig. 20(d).

Figure 19.

System identification model of Doppler waveform.

Figure 20.

Example of data sets for system identification.

4.4. System design of Vp waveform prediction

Figure 21 shows the valve-rejection algorithm using the coefficient sequences of Fig. 19 obtained by system identification. V_p and ECG (new data sets) were normalized and filtered. V_e (the prediction waveform of V_p) was generated by system prediction block. Based on differences between NVP and V_e, blending times of V_p and V_e were calculated. Based on blending times, the blending weights were changed. Blending weight generator was controlled so that NV_p becomes predominant in systolic phase. After blending process, V_m (the ideal waveform of V_p) was predicted.

Figure 22 shows the waveforms of Fig. 21. V_p and ECG are shown in Fig. 22(a) and (b), respectively. A complex heartbeat cycle waveform based on ECG is shown in Fig. 22(c). This indicates R-trigger timing, the systolic phase, and the diastolic phase by their amplitude level. V_p and V_e are shown in Fig. 22(d). Diff (absolute differences between V_p and V_e) waveform is shown in Fig. 22(e). Weighing functions (V_p-weight and V_e-weight) were generated based on Diff and the complex heartbeat cycle waveform. V_p-weight and V_e-weight, and V_m are shown in Fig. 22(f).

Figure 21.

Ideal V_m waveform prediction block diagram.

Figure 22.

Signal-processing waveforms of system prediction in Fig. 21.

4.5. Mathematical models

Parametric models such as ARMAX were used as mathematical models [10]. A parametric model is shown in formula (10).

$$P\left(z\right)\cdot y\left(z\right)=\frac{Q\left(z\right)}{R\left(z\right)}\cdot u\left(z\right)+\frac{S\left(z\right)}{T\left(z\right)}\cdot w\left(z\right)$$E10

Here, y(z) is an output, and u(z) and w(z) are an input and white noise, respectively. P(z), Q(z), R(z), S(z), and T(z) are the coefficient sequences. In Fig. 19, V_p and ECG correspond to y(z) and u(z), respectively. Several models that had different structures and orders were investigated.

4.6. Input data and mathematical model evaluation

Many clinical data sets of left ventricular outflow using a Doppler ultrasound system were acquired. Because outflow varies with individual differences, combined data from numerous volunteers was used for evaluation. Combined data has different waveforms, heartbeat cycles, and blood-flow sensitivities. Figure 23 shows combined V_p and ECG. Both V_p and ECG were sampled at 120 Hz sampling rate. The combined 16 heartbeat waveforms were used for simulations. Left ventricular outflows of volunteers A, B, and C are shown in Fig. 24(a)-(c) , respectively.

Several mathematical models were applied and evaluated. Orders of several models such as ARX model, ARMAX model, output error model (OE model), etc., were optimized, respectively. Next, the model fitnesses were evaluated by root-mean-square (RMS) errors. OE model was chosen for the valve-rejection algorithm because of the smallest RMS error.

Figure 23.

Combined data sets for system identification.

Figure 24.

Volunteers’ blood-flow data.

4.7. Verification

Finally, OE model was chosen as the optimal one for left ventricular outflow. In order to verify its performance of valve-rejection algorithm, the additional data other than volunteer A, B and C was needed. A different volunteer's data (data D) is shown in Fig. 25. Simulation results of the valve-rejection algorithm are shown in Fig. 25(b) using data D. During the first, second, third, and fifth heartbeat cycles, V_m traced ideal outflow. The valve signals were automatically rejected. But the valve signal was not sufficiently rejected during the fourth heartbeat cycle. Although there remains improvement of robustness, high performance of valve rejection was confirmed.

Figure 25.

Verification of valve rejection algorithm.

5. Conclusion

Based on the mathematical model that combined an ECG and biosignals (ultrasound Doppler image parameters, etc.), the system identification method to heart's blood flows was applied. With combination of the image parameter and the ECG, the effectiveness of a new gap-filling algorithm was confirmed. Moreover, with combination of the Doppler blood-flow waveform and the ECG, noises in heart's blood-flow measurement, such as valve regurgitation, were removed, and reliable automatic measurement of left ventricle outflow was realized. System identification using such a statistical method will be an important component for automatic measurement and diagnosis.