Ultrasound in Med. & Biol., Vol. 47, No. 4, pp. 10771088, 2021
Copyright © 2020 World Federation for Ultrasound in Medicine & Biology. All rights reserved.
Printed in the USA. All rights reserved.
0301-5629/$ - see front matter

https://doi.org/10.1016/j.ultrasmedbio.2020.12.011

 Original Contribution
MULTIFREQUENCY PHASED TRACKING METHOD FOR ESTIMATING VELOCITY
IN HEART WALL
TAGEDPYU OBARA,* SHOHEI MORI,y MOTOTAKA ARAKAWA,*,y and HIROSHI KANAI*,yTAGEDEN* Graduate School of Biomedical Engineering, Tohoku University, Sendai, Japan; and y Graduate School of Engineering, Tohoku
University, Sendai, Japan(Received 17 June 2020; revised 1 December 2020; in final from 12 December 2020)

Abstract—Local high-accuracy velocity estimation is important for the ultrasound-based evaluation of regional
myocardial function. The ultrasound phase difference at the center frequency of the transmitted signal has been
conventionally used for velocity estimation. In the conventional method, spatial averaging is necessary owing to
the frequency-dependent attenuation and interference of backscattered waves. Here, we propose a method forsuppressing these effects using multifrequency phase differences. The resulting improvement in velocity estimation in the heart wall was validated by in vivo experiments. In the conventional method, the velocity waveform
exhibits spike-like changes. The velocity waveform estimated using the proposed method did not exhibit suchchanges. Because the proposed method estimates myocardium velocity without spatial averaging, it can be used
for measuring heart wall dynamics involving thickness changes. (E-mail: mori@ecei.tohoku.ac.jp) © 2020
World Federation for Ultrasound in Medicine & Biology. All rights reserved.
Key Words: Ultrasound, Heart wall, Measurement of displacement, Phase difference, Interference, Thickness
change, Strain rate.

Guth et al. 1984). However, there are qualitative diagnosis methods that may be examiner dependent(Picano et al. 1991). Thus, it is necessary to develop
quantitative diagnosis methods based on the physiologic
mechanisms of myocardial contraction.
Many physiologic studies of myocardial contraction
have confirmed that the myocardium contracts owing toelectricalexcitation(Durrer
et
al.
1970;
McVeigh et al. 2002; Faris et al. 2003;
Ramanathan et al. 2004; Pernot et al. 2007;
Konofagou and Provost 2012; Grondin et al. 2019).
Therefore, the minute contraction response accompaniedby the conduction of electrical excitation has been measured in the time phase around the R-wave in the electrocardiogram, for regional myocardial function evaluation(Kanai and Tanaka 2011). When the myocardium contracts, its thickness increases owing to the volume conservation constraint. In previous studies in animalmodels using invasive methods, the regional myocardial
function was evaluated quantitatively by measuringchanges in myocardium thickness (Theroux et al. 1974;
Heyndrickx et al. 1975; Sasayama et al. 1976;
Sabbah et al. 1981; Hartley et al. 1983, 1991;
Homans et al. 1985). Therefore, it is possible to measure

INTRODUCTION
Ischemic heart diseases, such as myocardial infarction
and cardiac angina, are the leading factors determining
the number of years of life lost in 2016
(Foreman et al. 2018). Ischemic heart diseases negatively affect regional myocardial function. Therefore,
there is a need to develop highly accurate and repeatable
methods for evaluation of regional myocardial function;these diagnostic methods can help to prevent the onset ofischemic heart diseases and can increase the efficiencyof post-onset treatment. Ultrasound-based diagnostic
apparatuses are less expensive than other diagnostic
apparatuses, such as computed tomography and/or magnetic resonance imaging apparatuses, and can be used
for non-invasive real-time observation of organ dynamics. Given this, ultrasound-based medical evaluation
methods have been popular.
Conventionally, B-mode and M-mode imaging
have been used for evaluating myocardial contractile
function (Mason et al. 1979; Shapiro et al. 1981;
Address correspondence to: Shohei Mori, Graduate School ofEngineering, Tohoku University, Sendai 980-8579, Japan. E-mail:mori@ecei.tohoku.ac.jp

1077

1078

Ultrasound in Medicine & Biology

myocardial contraction by observing changes in myocardium thickness. Measurements of local changes in myocardium thickness (Kanai et al. 2006), myocardial strain(Sutherland
et
al.
2004;
Marwick
2006;
Grondin et al. 2015) and myocardial strain rate
(Kanai et al. 1997; Sutherland et al. 2004; Marwick 2006;
Yoshiara et al. 2007; Tanaka et al. 2014a, 2014b, 2017,2019) have been utilized for quantitative evaluation of
myocardial contractile function. Moreover, the strain
rate of the myocardium after the onset of ischemic heartdiseases has been studied, and the difference between
normal and infracted myocardia has been discussed(Jamal et al. 2002; Pislaru et al. 2004).
The changes in myocardium thickness, myocardial
strain and myocardial strain rate are calculated from velocities at multiple positions in the myocardium. In addition tomeasuring the myocardium contraction based on changes
in myocardium thickness, the minute contraction propagation velocity has been calculated from the delay times of
myocardial velocity waveforms (Matsuno et al. 2017;
Hayashi et al. 2019). Therefore, it is important to estimate
myocardial local velocity with high accuracy.
We have estimated the velocity in the ultrasoundbeam direction using the phased-tracking method(Kanai et al. 1996; Kanai et al. 1997). In the conventional
method, the quadrature detection method was applied to thereceived signal at the center frequency of the transmitted signal (Kanai et al. 1996, 1997, 2006; Yoshiara et al. 2007;
Tanaka et al. 2014a, 2014b, 2017; Matsuno et al. 2017;
Hayashi et al. 2019), and the velocity was estimated using
the phase difference between the received signals at the center frequency of the transmitted signal. However, the peakfrequency of the received signal does not correspond to thecenter frequency of the transmitted signal, as it depends on
two factors. One of these factors is the attenuation that
occurs when the ultrasound wave propagates through living
tissue. Owing to the frequency dependence of this attenuation, the signal’s high-frequency components are attenuatedmore than its low-frequency components especially for signals reflected from deep-tissue regions. Therefore, the peakfrequency of the received signal is lower than the center frequency of the transmitted signal. The other important factor
is the interference of backscattered waves from numerous
scatterers in the myocardium. Because this interference
depends on the scatterers’ spacing, the signal-to-noise ratio(SNR) may decrease at a certain frequency of the received
signal affected by the destructive interference. Thus, thepeak frequency of the received signal can be interference
dependent (Hasegawa and Kanai 2008a). Therefore, the conventional method is sensitive with respect to destructive
interference, which decreases the SNR at the center frequency of the transmitted signal.
Here, we propose a method for local and highly accurate estimation of myocardial velocity waveforms using theVolume 47, Number 4, 2021

phased-tracking method extended to multiple frequencies.
By calculating such multifrequency phase differences, evenwhen the SNR at a certain frequency is low, the phase differences at other frequencies (with high SNRs) can be used
to estimate velocity waveforms. A velocity estimator thatuses multifrequency phase differences has been consideredin several previous studies. One of them is the phase-sensitive 2-D motion estimator proposed by Hasegawa (2016).
The primary objective of that study was to expand a 2-D
motion estimator using the phase at a specific frequency(Salles et al. 2015) to one that uses multifrequency phases,
to improve the estimation accuracy of the displacement in
the lateral direction, that is, perpendicular to the beam axis,
using a smaller spatial window. However, the effects ofattenuation and interference of scattered waves have notbeen studied thoroughly. In the present study, the main purpose of the proposed method is to suppress the influence of
the attenuation and interference of scattered waves in theheart wall dynamics measurements, where myocardium
thickness changes owing to myocardium contraction.
We compared the proposed phased-tracking method
extended to multiple frequencies with the conventional
phased-tracking method that uses the phase difference at
the center frequency of the transmitted signal, and evaluated the reproducibility of the velocity waveforms estimated using the proposed method. We discuss theusefulness of the proposed method for measuring heart
wall dynamics involving thickness changes.

METHODS
Conventional methodThe
conventional
phased-tracking
method(Kanai et al. 1996, 1997) at the center frequency of thetransmitted signal is explained below and in Figure 1. (I)Fig. 1. Schematic of the estimation of velocity using the phasegradient of the cross-spectrum in the conventional method(blue lines) and in the proposed method (red lines).
RF = radiofrequency.

Estimating velocity in heart wall  Y. OBARA et al.

The received radiofrequency (RF) signal sn ðtÞ aroundthe analysis point in the nth frame is cut out using a window function. (II) The complex frequency spectrum Sn ðf Þ
is obtained using the discrete Fourier transform (DFT).
The DFT is also applied to the received signal snþ1 ðtÞ ofthe ðn þ 1Þth frame. (III) The cross-spectrum Cn ðf Þ
between the consecutive frames is calculated as
Cn ðf Þ ¼ Sn ðf Þ ¢ Snþ1 ðf Þ

ð1Þ

where  denotes the complex conjugate. The amplitudeof the cross-spectrum, jCn ðf Þj, represents the correlation,in the frequency domain, between the frames of thereceived signal, and the phase ﬀCn ðf Þ of the cross-spectrum represents the phase difference between the frames^ d ðnÞ isof the received signal. The estimated velocity v
calculated as^ d ðnÞ ¼v

c0 fFR ﬀCn ðfc Þ
¢
fc
4p

ð2Þ

where c0 is the speed of sound in the living tissue, fFR is
the frame rate and fc is the center frequency of the transmitted signal.

Proposed method
The proposed phased-tracking method, which isextended to multiple frequencies, is described below.
Assuming that the velocity vd ðnÞ between the frames is
uniform in the window of the DFT and the phase ﬀCn ðf Þof the cross-spectrum is linear, the phase gradient an of
the cross-spectrum is expressed as
an ¼ 

4p¢ vd ðn Þ
c0 fFR

ð3Þ

Therefore, the velocity vd ðnÞ between the consecutive
frames can be calculated from the phase gradient an of
the cross-spectrum.
In the proposed method, the phase gradient an of thecross-spectrum is estimated from the phase ﬀCn ðf Þ withmultiple frequencies. Assuming that the cross-spectrum
phase is linear, the root mean squared error (RMSE) of theestimated phase gradient, an ðaÞ, weighted by the crossspectrum amplitude jCn ðf Þj, is defined asvﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
uPfs
2
u 2
u f ¼ 0 jCn ðf Þj ¢ jﬀCn ðf Þ  af jtan ðaÞ ¼
ð4Þ
Pf2s
jCn ðf Þj
f ¼ 0
for the sampling frequency fs and arbitrary variable a.
When the value of af is greater than (less than) 180˚, 360˚is subtracted (added) from (to) the value of af to deal with^ d ðnÞ can be estimatedaliasing. From eqn (3), the velocity v^ n of the cross-spectrum that minifrom the phase gradient amizes the RMSE an ðaÞ. The velocity is calculated as^ d ðnÞ ¼v

c0 fFR
^nÞ
¢ ða
4p

^ n ¼ arg min an ðaÞ
a
a

1079

ð5Þ
ð6Þ

The concept underlying this estimation method is the sameas that underlying the phase-sensitive 2-D motion estimator
proposed by Hasegawa (2016); however, the main purpose
is different, as mentioned in the Introduction.
Figure 1 schematically illustrates estimation of thevelocity using the phase gradient of the cross-spectrum inthe conventional and proposed methods. The blue lines represent the spectra and phases for the conventional method,whereas the red lines represent the spectra and phases forthe proposed method. In the conventional method, the estimation accuracy of the cross-spectrum phase might belower, owing to the influence of attenuation and the interference, because only the phase at the center frequency of
the transmitted signal is used. Because the higher-frequency components of the ultrasound wave are more attenuated, the peak frequency fp of the signal received from theheart wall is lower than the center frequency fc of the transmitted signal, as illustrated in Figure 1. Therefore, in theconventional method, the velocity might be estimated fromthe phase difference at the frequency with a low SNR of
the cross-spectrum.
Furthermore, the phase difference at a certain frequency may be affected by the interference, depending on
the position of the myocardium to be analyzed. Because
there are numerous scatterers in the myocardium that aresmaller than the resolution of the ultrasound wave, thereceived signal is the interference of numerous backscatteredwaves. When destructive interference occurs, for example,
as a result of the opposite phase in assuming two scatterers,
the SNR decreases dramatically, and the cross-spectrum
phase cannot be calculated accurately. When destructive
interference occurs at the center frequency of the transmitted
signal, the estimation accuracy of the velocity using the conventional method with the phase difference at the center frequency of the transmitted signal becomes lower. In the^ n isproposed method, the cross-spectrum phase gradient aestimated by using the cross-spectrum amplitude jCn ðf Þj as
the weighting function. Therefore, the influence of the phasedifference at the frequency at which the SNR is lower owing
to the attenuation and/or interference can be suppressed.
Measurement environment
In vivo measurements were made on the interventricular septum (IVS) of healthy subjects in their early twenties,
using ultrasound diagnostic apparatuses. The measurement
was approved by the Ethics Committee of Tohoku University, and the subjects agreed to participate in this study. For
measurement of subjects A, B and C, a plane wave with the1080

Ultrasound in Medicine & Biology

center frequency fc of 3.75 MHz was transmitted from a sector probe (UST-52101 N; Hitachi-Aloka-Medical Ltd,Tokyo, Japan) connected to an ultrasonic diagnostic apparatus (Prosound a-10; Hitachi-Aloka-Medical Ltd, Tokyo,
Japan), and high-frame-rate measurements (fFR = 860 Hz)
were made using parallel beamforming (Hasegawa andKanai 2008b). The sampling frequency fs was 15 MHz, andthe sound velocity c0 in the living tissue was assumed to be
1540 m/s; thus, the depth interval of the sampling points was
0.051 mm. For measurement of subjects D, E and F, afocused wave with the center frequency fc of 3.75 MHz wastransmitted from a sector probe (UST-52101; HitachiAloka-Medical Ltd, Tokyo, Japan) connected to an ultrasonic diagnostic apparatus (SSD6500; Hitachi-Aloka-Medical Ltd, Tokyo, Japan). By reducing the number of beams,
high-frame-rate measurement (fFR  630 Hz) was realized.
The sampling frequency fs was 20 MHz, and the soundvelocity c0 in the living tissue was assumed to be 1540 m/s;
thus, the depth interval of the sampling points was
0.038 mm.
In these measurements, the electrocardiogram(ECG) and phonocardiogram (PCG) were acquired
simultaneously. For each subject, the ECG was measured using the three-point lead method, and the ECG
waveform was obtained for lead II. The PCG waveform
was measured using a small microphone attached to thesubject’s chest.
Figure 2 is the initial-frame B-mode image of a subject and the analysis point where the velocity was estimated (yellow point). Around the time phase of the ECG
R-wave representing the onset of the minute contraction
response, the velocity at the analysis point was estimatedusing the phased-tracking method at the center frequencyof the transmitted signal (conventional method, eqn ð2Þ)and using multiple frequencies (proposed method, eqn
ð5Þ). In the velocity estimation, the received RF signal

Fig. 2. B-mode image of the left ventricle in the long-axis viewand the analyzed point at which the velocity was estimated (yellow
point). ECG = electrocardiogram; PCG = phonocardiogram.

Volume 47, Number 4, 2021

was cut out using the Hanning window in the range of§1.22 mm around the analysis point in the depth direction. Then, the cross-spectrum Cn ðf Þ in eqn (1) was
obtained, and the velocity between two consecutive
frames was estimated. The analysis point was moved
with the movement of IVS. The movement of IVS was
measured by integrating the estimated velocity waveform.
To compare the proposed method with the conventional method using the RMSE an ðaÞ in eqn (4), there
was a need to calculate the RMSE in the conventional
method. Therefore, the phase gradient of the cross-spec^ n in eqn (4) was taken place by ﬀCn ðfc Þ=fc to caltrum aculate the RMSE for the conventional method using thephase at the center frequency of the transmitted signal.

RESULTS
Figure 3 illustrates (a) the ECG waveform, (b) the
PCG waveform, (c) the velocity waveform and (d) the
RMSE in the phase gradient estimation. On the verticalaxis in Figure 3c, a downward positive value indicates thecontraction movement toward the left ventricle (LV), andan upward negative value indicates the movement toward
the right ventricle (RV). The dotted lines represent the
maximal velocities that can be measured using the phasedifference at the center frequency fc of the transmitted signal. The horizontal axis is the time from the ECG R-wave.
Figure 3c, 3d illustrates the results obtained using thephased-tracking method extended to multiple frequencies(proposed method, red lines) and using the method thatuses only the center frequency fc of the transmitted signal(conventional method, blue lines).

Fig. 3. (a) Electrocardiogram (ECG) waveform, (b) phonocardiogram (PCG) waveform, (c) velocity waveform and (d) minimal root mean square error (RMSE) of the phase gradientestimation, where blue lines represent the results obtainedusing the conventional method, which uses the phase difference at the center frequency of the transmitted signal, and redlines represent the results obtained using the proposed methodthat uses phase differences with multiple frequencies. LV = left
ventricle; RV = right ventricle.

Estimating velocity in heart wall  Y. OBARA et al.

The velocity waveform estimated using the conventional method exhibited spike-like changes, while the waveform estimated using the proposed method did not exhibit
spike-like changes. Moreover, when the velocity waveform
estimated using the conventional method exhibited spikelike changes, the RMSE in the conventional method was
extremely large compared with that in the proposed method.
This suggests that the spike-like velocity changes observed
for the estimation using the conventional method could be
attributed to the erroneous estimation.
Let us discuss the cause of the spike-like velocitychange in the velocity waveform estimated using theconventional method. We focused on the detail of thevelocity estimation for the n1 th frame (44.9 ms from the
ECG R-wave) with the maximal RMSE for the conventional method in Figure 3d, at which the positive spikelike velocity change was estimated as illustrated in
Figure 3c. Figure 4 illustrates (a) the received RF signals
before cutting out using the Hanning window, (b) the
amplitudes of the complex frequency spectra, jSn ðf Þj, (c)the phases of the complex frequency spectra, ﬀSn ðf Þ, (d)the amplitude of the cross-spectrum, jCn ðf Þj, and (e) thephase of the cross-spectrum, ﬀCn ðf Þ. In Figure 4ac, theblue lines correspond to the frequency spectrum Sn ðf Þ ofthe frame with the maximal RMSE, for the conventional
method, while the red lines correspond to the frequency
spectrum Snþ1 ðf Þ in the next frame. In Figure 4d, 4e, thesolid lines represent the cross-spectrum Cn ðf Þ between
the nth and ðn+1)th frames, the point (✕) represents the
cross-spectrum phase ﬀCn ðfc Þ at the center frequency fc
of the transmitted signal, for the conventional method,
and the red dashed line represents the gradient of the^ n estimated using the proposedcross-spectrum phase amethod with multiple frequencies.

1081

In Figure 4b, the peak frequency of the received signal is lower than the center frequency fc of the transmitted signal. This is attributed to the frequency-dependent
attenuation. Despite the frequency-dependent attenuation, jSn ðfc Þj and jSnþ1 ðfc Þj are lower than the amplitudeof the complex frequency spectrum at the adjacent frequency, which is approximately 4 MHz. This is attributed to the destructive interference. Although theanalyzed position in the myocardium was the same, jSn ðfÞj differs from jSnþ1 ðf Þj. This is attributed to the change
in the interference state. In Figure 4d, the peak frequency
fp of the cross-spectrum amplitude jCn ðf Þj is approximately 3 MHz. The cross-spectrum amplitude at the center frequency fc of the transmitted signal jCn ðfc Þj is37 dB in Figure 4d, relative to that at the peak frequency fp of the received signal. In Figure 4e, the phaseof the cross-spectrum, ﬀCn ðfc Þ, at the center frequency fc
of the transmitted signal is significantly different fromthe phase at the peak frequency. These results indicatethat the SNR decreased at the center frequency fc of thetransmitted signal owing to the influence of the attenuation and the interference, and ﬀCn ðfc Þ used in the conventional method was not calculated as the phase difference
reflecting the actual heart wall motion. On the other
hand, the attenuation and interference effects were suppressed in the proposed method, because the estimated^ n was calculated from the phasecross-phase gradient a
difference with multiple frequencies using the crossspectrum amplitude jCn ðf Þj as a weighting function. In
Figure 4d, the 20-dB bandwidth in the velocity estimation was 2.453.37 MHz.
We evaluated the reproducibility of the velocity
waveforms estimated using the proposed method. Figures 57 illustrate the results around the time phase of

Fig. 4. (a) Received radiofrequency signal sn1 ðtÞ, (b) amplitude of the complex frequency spectrum, jSn1 ðf Þj, and (c)phase of the complex frequency spectrum, ﬀSn1 ðf Þ, where the blue line represents Sn1 ðf Þ; where RMSE was maximal forthe conventional method, and the red line represents Sn1 þ1 ðf Þ. (d) The amplitude of the cross-spectrum, jCn1 ðf Þj, of then1 th frame with maximal root mean square error (RMSE) in the conventional method and (e) phase of the cross-spec^ n1 , estimated using the protrum, ﬀCn1 ðf Þ, where the dashed line represents the phase gradient of the cross-spectrum, a
posed method, and the symbol ✕ indicates the cross-spectrum phase ﬀCn1 ðfc Þ at the center frequency of the transmitted
signal, using the conventional method. ECG = electrocardiogram; PCG = phonocardiogram.

1082

Ultrasound in Medicine & Biology

Volume 47, Number 4, 2021

Fig. 5. Results of the estimated velocity waveform for subject A in two consecutive heartbeats. (a) Electrocardiogram (ECG)
waveform. (b) Phonocardiogram (PCG) waveform. (c) Velocity waveform estimated using the conventional method that uses
the phase difference at a single center frequency of the transmitted signal. (d) Velocity waveform estimated using the proposed
method that uses phase differences with multiple frequencies. LV = left ventricle; RV = right ventricle.

Fig. 6. Results of the estimated velocity waveform for subject B in two consecutive heartbeats. (a) Electrocardiogram(ECG) waveform. (b) Phonocardiogram (PCG) waveform. (c) Velocity waveform estimated using the conventional
method that uses the phase difference at the center frequency of the transmitted signal. (d) Velocity waveform estimatedusing the proposed method that uses phase differences with multiple frequencies. LV = left ventricle; RV = rightventricle.

the ECG R-wave for two consecutive heartbeats and for
different subjects. The velocity waveforms estimatedusing the proposed method are provided in Figures 5d,
6d and 7d and are quite similar. On the other hand, thevelocity waveforms estimated using the conventional
method are provided in Figures 5c, 6c and 7c and arequite different because spike-like velocity changes
owing to the lower SNR occur independently from theheart wall dynamics. In Figure 7c, the velocity waveforms estimated using the conventional method exhibit
steep discontinuity, which is caused by aliasing after
20 ms from the R-wave. This is because the velocity inthe heart wall exceeded the maximal velocity that can bemeasured using the phase difference at the center

frequency fc of the transmitted signal. In the proposed
method, the velocity was estimated using the phase gradient of the cross-spectrum. In the proposed method, if
the change in the cross-spectrum phase between adjacentfrequency components was less than 180˚, the phase gradient of the cross-spectrum can be calculated without aliasing. In the present measurement, the change in the
cross-spectrum phase between the peak frequency andits adjacent frequency was approximately 10˚. Therefore,
in the proposed method, the velocity was measured without aliasing. For measurement of subjects D, E and F,
the reproducibility of the velocity waveforms was evaluated by calculating the standard deviation (SD) in four
consecutive heartbeats. The values of the SD averaged

Estimating velocity in heart wall  Y. OBARA et al.

1083

Fig. 7. Results of the estimated velocity waveform for subject C in two consecutive heartbeats. (a) Electrocardiogram(ECG) waveform. (b) Phonocardiogram (PCG) waveform. (c) Velocity waveform estimated using the conventional
method that uses the phase difference at the center frequency of the transmitted signal. (d) Velocity waveform estimatedusing the proposed method that uses phase differences with multiple frequencies. LV = left ventricle; RV = rightventricle.

over frames around the time phase of the ECG R-wave
were 7.3 (conventional) and 1.8 (proposed) in subject D,
7.7 (conventional) and 3.0 (proposed) in subject E, and
6.5 (conventional) and 2.8 (proposed) in subject F,
respectively. In three subjects, the SD in the proposed
method was smaller than that in the conventional
method. Therefore, we concluded that the velocity estimation using the proposed method was more precise
than that obtained by the conventional method. For subjects A, B and C, the SD was not calculated because the
RF data were acquired in only two heartbeats.
In the velocity waveforms estimated using the proposed method, the positive velocity toward the LV wasobserved during 30 to 0 ms from the ECG R-wave, as
illustrated in Figures 5d, 6d and 7d. In our previous studies, it was suggested that this velocity was related to theminute contraction response caused by electrical excitation (Kanai and Tanaka 2011; Matsuno et al. 2017;
Hayashi et al. 2019). In our future work, we will validate
the relationship between the velocity waveforms and theactual electrical excitation. The precise velocity estimation by the proposed method would contribute to revealing this relationship.
DISCUSSION
Improvement of velocity estimation using phasedifferences with multiple frequencies
The proposed method improved velocity estimation
because it accounted for the effects of attenuation andinterference, which are not considered in the conventional method.
In the conventional phased-tracking method thatuses the phase at the center frequency fc of thetransmitted signal, the estimation accuracy is negatively
affected by the frequency-dependent attenuation. In
Figure 4d, the peak frequency fp of the received signal islower than the center frequency fc of the transmitted signal. This is because higher-frequency components weremore affected by the attenuation when propagating in
the living tissue. Therefore, the SNR at the center frequency of the transmitted signal was lower, and thevelocity estimation accuracy was lower in the conventional method, which used only the phase at the center
frequency of the transmitted signal.
The effect of the interference with other backscattered waves should be discussed in the measurement ofheart wall measurements. The main scatterers in the
myocardium are collagen fibers. The myocardium has alayered structure, consisting of several myocardial fibers,
and there are abundant collagen fibers between thelayers. Because collagen fibers are smaller than the ultrasound spatial resolution, the signals received from the
myocardium with numerous scattered waves interferewith each other, and the B-mode image exhibits aspeckle pattern. The scattered waves can build up (constructive interference, or in-phase waves) or cancel out(destructive interference, or out-of-phase waves). Therefore, the specific frequency component of the cross-spectrum decreases depending on the spacing among thenumerous scatterers in the myocardium. In Figures 4band 4d, the frequency spectrum amplitude jSn ðf Þj, jSnþ1 ðf Þj and the cross-spectrum amplitude jCn ðf Þj
decreased at the center frequency of the transmitted signal. This occurred owing to the destructive interferenceat the analyzed position, especially at the center frequency of the transmitted signal. Thus, the phase difference calculated at the center frequency of the1084

Ultrasound in Medicine & Biology

transmitted signal would be affected by the backgroundnoise.
A notable disadvantage of the conventional method
is not only to use the phase difference at the center frequency of the transmitted signal but also to use the phasedifference for only a single frequency of the received
signal. In the heart wall, rapid changes in the wall thickness associated with myocardial contraction change thespacing among the scatterers. Because the amplitude of
the myocardial strain rate exceeds 4 (m/s)/m when measured in the left ventricular long-axis view
(Yoshiara et al. 2007; Tanaka et al. 2014a, 2014b, 2017),
an approximately 10-mm change in thickness in the spatial window occurs over 1 ms; thus, changes in the interference state between two consecutive frames occur
even for high-frame-rate measurements. Even if the
SNR is not low, a change in the interference between
two consecutive frames induces a phase difference thatis different from that caused by the actual heart wallmotion. Because of the state of the interference changes
owing to rapid changes in heart wall thickness, the conventional phased-tracking method (Kanai et al. 1996,1997) and the phased-tracking method that was developed for addressing arterial wall thickness changes(Hasegawa et al. 2002; Hasegawa and Kanai 2006)
would be difficult to apply to heart wall measurements.
In the proposed method, the phased-tracking
method extended to multiple frequencies alleviated thenegative effects of attenuation and interference on estimation accuracy. Because the velocity is estimated fromthe phases with multiple frequencies using the crossspectrum amplitude jCn ðf Þj as a weighting function, theinfluence of the phase at the frequency where the SNR
decreases owing to attenuation and interference isVolume 47, Number 4, 2021

alleviated. Furthermore, the influence of the change in
the interference state is alleviated when estimating thevelocity based on the phase differences with multiple
frequencies. If the phase difference caused by the change
in the interference state is random at different frequencies, combining the phase differences for several frequencies is likely to be useful for suppressing the effectof the change in the interference state. As a result,smooth waveforms (without spike-like velocity changes)
were obtained with high precision.
Figure 8 illustrates (a) the ECG waveform, (b) the
PCG waveform, (c) the velocity waveform, (d) the
RMSE of the phase gradient estimation and (e) the depthof the tracking point during a cardiac cycle for subject
D. In the conventional method, the depth of the tracking
point did not return to its original position after one cardiac cycle. Because the displacement is obtained by integrating the velocity waveform, an erroneous velocity
estimation leads to an unrepairable error in tracking theanalyzed point. The deviation from the original position
after one cardiac cycle was improved by the proposed
method compared with the conventional method in subject D, as illustrated in Figure 8e, and similar results
were confirmed in other five subjects. Therefore, we concluded that the proposed method has the potential to estimate the velocity with higher accuracy than theconventional method. In our future work, the accuracyof velocity estimation using the proposed method will be
evaluated by the simulation or the experiment using aphantom which has a local change in thickness.
Usefulness in measurements of myocardial contraction
We examined the usefulness of the proposed method
for estimating the velocity for measuring the heart wall’s

Fig. 8. Results of the estimated velocity waveform for subject D during a cardiac cycle. (a) Electrocardiogram (ECG)
waveform. (b) Phonocardiogram (PCG) waveform. (c) Velocity waveform. (d) Minimal root mean square error (RMSE)of the phase gradient estimation. (e) Depth of the analyzed point, where blue lines represent the results obtained using
the conventional method that uses the phase difference at the center frequency of the transmitted signal, while the redlines represent the results obtained using the proposed method that uses phase differences with multiple frequencies.

Estimating velocity in heart wall  Y. OBARA et al.

contraction response. The myocardial contraction response
was measured by calculating the myocardial strain rate
(Kanai et al. 1997; Sutherland et al. 2004;
Yoshiara et al. 2007; Tanaka et al. 2014a, 2014b, 2017)
and the minute contraction propagation speed from thevelocity waveforms at multiple heart wall locations(Matsuno et al. 2017; Hayashi et al. 2019). However,
because the velocity waveforms obtained using the conventional method are affected by attenuation and interference,as mentioned above, it was difficult to reliably estimate the
myocardial strain rate and the minute contraction propagation velocity from the velocity estimated at a single measurement point. For this reason, when measuring thecontraction response using the conventional method, thephase ﬀCn ðfc Þ at the center frequency of the transmitted signal was spatially averaged around the analysis point in thedepth direction, and the cross-correlation function was used
for estimating the spatially averaged velocity(Kanai et al. 1996, 1997).
The usefulness of the proposed method in measuring
the change in thickness of the heart wall was validated bymeasuring the transmural thickness of the heart wall using
the conventional and proposed methods during a cardiac
cycle in six healthy subjects. The change in wall thickness
was calculated from the differences between the velocity
waveforms at the LV side and the RV side of the IVS.
Figure 9 illustrates results for subject D with (a) the ECG
waveform, (b) the PCG waveform, (c) the velocity waveforms estimated using the conventional method with spatial
averaging, (d) the velocity waveforms estimated using the
proposed method and (e) the thicknesses of the heart wall
for subject D. The results in Figure 9c, 9d illustrate the positive velocity waveforms toward the LV (contraction)

1085

during the ejection period, and then the negative velocity
waveforms toward the RV during the filling period.
Because the amplitude of the velocity at the LV side waslarger than that at the RV side during the ejection and filling periods, the thickness increased during the ejection
period and decreased during the filling period, as illustrated
in Figure 9e. Around the time phase of the ECG R-wave,
there was almost no change in the thickness.
The waveform of the thickness measured by the
proposed method was similar to that measured using
ultrasonic
crystalsintheprevious
studies(Sasayama et al. 1976; Hartley et al. 1983, 1991). The
thickness measured using the proposed method (red line)
returned to the original thickness after one cardiac cycle,
but that measured using the conventional method (blackline) did not. These results were confirmed in six subjects. In the six healthy subjects, the thickening fraction(Hartley et al. 1983) was evaluated. The thickening fractions using the proposed method were in the range
20%38%. The values of thickening fraction agreed
with those in previous studies using ultrasonic crystals:
30% in dogs (Sasayama et al. 1976), 23%38% in piglets (Hartley et al. 1983) and 20% in dogs
(Hartley et al. 1991), respectively. Therefore, it was suggested that the proposed method was useful in measuring
the change in the thickness of the heart wall.
The averaging effect of the velocity is discussed in
Figure 10. The black line in Figure 10b represents thevelocity waveform obtained using the conventional methodwith spatial averaging. In the conventional method withspatial averaging (black line in Fig. 10b), the cross-spectraat the center frequency of the transmitted signal were averaged in the range of §1.22 mm from the analysis point in

Fig. 9. Results of the estimated velocity waveform and the thickness of the heart wall for subject D during a cardiac
cycle. (a) Electrocardiogram (ECG) waveform. (b) Phonocardiogram (PCG) waveform. (c) Velocity waveform estimatedusing the conventional method that uses the phase difference at a single center frequency of the transmitted signal withspatial averaging for cross-correlation function. (d) Velocity waveform estimated using the proposed method that uses
phase differences with multiple frequencies. (e) Thickness of heart wall estimated using the distance from the depth of
the analyzed point at the right ventricle (RV) side to that at the left ventricle (LV) side.

1086

Ultrasound in Medicine & Biology

Volume 47, Number 4, 2021

Fig. 10. (a) Electrocardiogram (ECG) waveform. (b) Velocity waveform. The blue line represents the results obtainedusing the conventional method that uses the phase difference at the center frequency of the transmitted signal. The blackline represents the averaged results for velocities in the depth direction obtained using the conventional method that uses
the phase difference at the center frequency of the transmitted signal. The red line represents the results obtained using
the proposed method that uses phase differences with multiple frequencies. (c) Conceptual diagram of the calculation of
the myocardial strain rate using the conventional method. (d) Conceptual diagram of the calculation of the myocardial
strain rate using the proposed method. Yellow circles represents the analyzed points at which the velocity was estimated.
DFT = discrete Fourier transform; LV = left ventricle; RV = right ventricle.

the depth direction. As illustrated in Figure 10b, spike-likevelocity changes estimated using the conventional method
were alleviated by spatial averaging.
There are differences between the velocity waveform estimated using the conventional method with spatial averaging and that estimated using the proposed
method without spatial averaging. This is attributed tochanges in the heart wall thickness that are associatedwith myocardial contractions. If the velocity is uniform,
that is, if the thickness does not change, within the rangeof averaging, the accuracy of the velocity estimation at
the analysis point is improved by averaging. However,
because the heart wall thickness changes locally owing
to contractions, the velocities in the averaging range inthe depth direction differ depending on the depth, andthe accuracy of the velocity estimation at the analysis
point decreases with averaging. Because the proposed
method can yield smooth velocity waveforms without
averaging, the velocity can be estimated locally, forexample, at a single point, compared with the conventional method. Thus, the proposed velocity estimation
method with multiple frequencies is useful for measuring
heart wall dynamics that involve thickness changes.
Finally, let us discuss how the proposed method
improves the locality in measuring the heart wall dynamics by considering the myocardial strain rate measurement as an example. The myocardial strain rate of thenth frame, SRðnÞ, is calculated using the velocity vdþ ðnÞat the deep analysis point, the velocity vd ðnÞ at the shallow analysis point and the distance h0 between their two
points (Yoshiara et al. 2007; Tanaka et al. 2014a, 2014b,
2017), asSR ðnÞ ¼

vdþ ðnÞ  vd ðnÞ
:
h0

ð7Þ

In previous studies, the distance h0 between the twopoints was only 821 mm. Figure 10c, 10d illustrates theconceptual diagram of ranges of velocities used for calculating the myocardial strain rate in the conventional andproposed methods. Because the analyzed region in the conventional method was averaged in the depth direction, thevelocity estimated at the analysis point included different
velocities in the averaging window, as illustrated in
Figure 10c. When the DFT window is §1.22 mm for thevelocity estimation, the overlaps between the two windows
in the myocardial strain rate measurements are 83% and
66% for the conventional and proposed methods, respectively. The proposed method estimated the velocity at theanalysis point in a more local manner. Therefore, the myocardial strain rate is calculated more locally in the proposed
method, compared with the conventional method. In our
future work, the local strain rate will be measured using the
proposed method.

Estimating velocity in heart wall  Y. OBARA et al.

Clinical application
The proposed method is applicable to existing ultrasonic diagnostic apparatuses with only the software
upgrade for use offline. For use in real time, a hardware
board may be a need if the computational capacity is insufficient in the existing ultrasonic diagnostic apparatus. Thisvelocity estimation procedure can be achieved by digital
signal processing so it is not difficult to implement our proposed method in existing ultrasonic diagnostic apparatuses.
In the clinical application, first, the region of interest (ROI)
is set on the heart wall, and the velocities in the top and thebottom of the ROI are estimated using the proposed
method. Then, the thickening (strain) or strain rate in the
ROI is calculated from the velocities.
The velocity estimated using the proposed method has
the parallel component to the ultrasound beam. In the previous study, it was confirmed that the strain and strain rate
calculated from the velocities in the ultrasound beam direction have the potential to quantify myocardial ischemia andto evaluate the LV mass regression after treatments(Sutherland et al. 2004). It is possible to expand the proposed method to the 2-D motion estimator using the 2-D
Fourier transform (Hasegawa 2016), if needed.
CONCLUSIONSIn the present study, we proposed a velocity estimation method that alleviates the effects of attenuation andinterference, which negatively affect the velocity estimation in the conventional method that uses the phase difference at the center frequency of the transmitted signal.
Using the proposed method, in which the velocity wasestimated using the phase differences with multiple frequencies, we were able to obtain the velocity waveformin the direction of the ultrasound beam in the heart wallwith high precision. Moreover, it was confirmed that theestimation precision of the conventional velocity estimation method using the phase difference at the center frequency of the transmitted signal decreased owing to thefrequency-dependent attenuation and interference of
scattered waves.
The proposed method can estimate the local velocity waveform of the myocardium without averaging in
the depth direction, while the conventional method
requires averaging. This locality of velocity estimation
is essential for measuring the dynamics of the heart wallthat exhibits thickness changes because the velocity
depends on the analyzed point.
Acknowledgments—This work was supported in part by Japan Society
for the Promotion of Science (JSPS) KAKENHI 19 K22943.
Conflict of interest disclosure—The authors declare that they have nocompeting interests.

1087

REFERENCES
Durrer D, van Dam RT, Freud GE, Janse MJ, Meijler FL, Arzbaecher
RC. Total excitation of the isolated human heart. Circulation
1970;41:899–912.
Faris OP, Evans FJ, Ennis DB, Helm PA, Taylor JL, Chesnick AS, Guttman MA, Ozturk C, McVeigh ER. Novel technique for cardiac electromechanical mapping with magnetic resonance imaging tagging andan epicardial electrode sock. Ann Biol Eng 2003;31:430–440.
Foreman KJ, Marquez N, Dolgert A, Fukutaki K, Fullman N,
McGaughey M, Pletcher MA, Smith AE, Tang K, Yuan CW,
Brown JC, Friedman J, He J, Heuton KR, Holmberg M, Patel DJ,
Reidy P, Carter A, Cercy K, Chapin A, Douwes-Schultz D, FrankT, Goettsch F, Liu PY, Nandakumar V, Reitsma MB, Reuter V,
Sadat N, Sorensen RJD, Srinivasan V, Updike RL, York H, Lopez
AD, Lozano R, Lim SS, Mokdad AH, Vollset SE, Murray CJL.
Forecasting life expectancy, years of life lost, and all-cause andcause-specific mortality for 250 causes of death: Reference andalternative scenarios for 201640 for 195 countries and territories.
Lancet 2018;392:2052–2090.
Grondin J, Wan E, Gambhir A, Garan H, Konofagou EE. Intracardiac
myocardial elastography in canines and humans in vivo. IEEE
Trans Ultrason Ferroelectr Freq Control 2015;62:337–349.
Grondin J, Wang D, Grubb CS, Trayanova N, Konofagou EE. 4D cardiac electromechanical activation imaging. Comput Biol Med
2019;113:1–14.
Guth B, Savage R, White F, Hagan A, Samtoy L, Bloor C. Detection ofischemic wall dysfunction: Comparison between M-mode echocardiography and sonomicrometry. Am Heart J 1984;107:449–457.
Hartley CJ, Latson LA, Michael LH, Seidel CL, Lewis RM, Entman
ML. Doppler measurement of myocardial thickening with a single
epicardial transducer. Am J Physiol 1983;245:HlO66–H1072.
Hartley CJ, Litowitz H, Rabinovitz RS, Zhu W-X, Chelly JE, Michael
LH, Bolli R. An ultrasonic method for measuring tissue displacement: Technical details and validation for measuring myocardial
thickening. IEEE Trans Biomed Eng 1991;38:735–747.
Hasegawa H. Phase-sensitive 2D motion estimators using frequency
spectra of ultrasonic echoes. Appl Sci 2016;6:195.
Hasegawa H, Kanai H. Improving accuracy in estimation of artery-walldisplacement by referring to center frequency of RF echo. IEEE
Trans Ultrason Ferroelectr Freq Control 2006;53:52–63.
Hasegawa H, Kanai H. Reduction of influence of variation in center
frequencies of RF echoes on estimation of artery-wall strain. IEEE
Trans Ultrason Ferroelectr Freq Control 2008a;55:1921–1934.
Hasegawa H, Kanai H. Simultaneous imaging of artery-wall strain andblood flow by high frame rate acquisition of RF signals. IEEE
Trans Ultrason Ferroelectr Freq Control 2008b;55:2626–2639.
Hasegawa H, Kanai H, Koiwa Y. Modified phased tracking method formeasurement of change in thickness of arterial wall. Jpn J Appl
Phys 2002;41:3563–3571.
Hayashi A, Mori S, Arakawa M, Kanai H. Local Two-dimensional distribution of propagation speed of myocardial contraction for ultrasonic visualization of contraction propagation path. Jpn J Appl
Phys 2019;58:SGGE05.
Heyndrickx GR, Millard RW, McRitchie RJ, Maroko PR, Vatner SF.
Regional myocardial functional and electrophysiological alterations after brief coronary artery occlusion in conscious dogs. J Clin
Invest 1975;56:978–985.
Homans DC, Asinger R, Elsperger KJ, Erlien D, Sublett E, Mikell F,
Bache RJ. Regional function and perfusion at the lateral border ofischemic myocardium. Circulation 1985;71:1038–1047.
Jamal F, Kukulski T, Sutherland GR, Weidemann F, D’hooge J, Bijnens B, Derumeaux G. Can changes in systolic longitudinal deformation quantify regional myocardial function after an acuteinfarction? An ultrasonic strain rate and strain study. J Am Soc
Echocardiogr 2002;15:723–730.
Kanai H, Tanaka M. Minute mechanical-excitation wave-front propagation in human myocardial tissue. Jpn J Appl Phys 2011;50 07HA01.
Kanai H, Sato M, Koiwa Y, Chubachi N. Transcutaneous measurement
and spectrum analysis of heart wall vibrations. IEEE Trans Ultrason Ferroelectr Freq Control 1996;43:791–810.

1088

Ultrasound in Medicine & Biology

Kanai H, Hasegawa H, Chubachi N, Koiwa Y, Tanaka M. Noninvasive
evaluation of local myocardial thickening and its color-coded
imaging. IEEE Trans Ultrason Ferroelectr Freq Control
1997;44:752–768.
Kanai H, Hasegawa H, Imamura K. Spatial distribution measurementof heart wall vibrations generated by remote perturbation of innerpressure. Jpn J Appl Phys 2006;45:4718–4721.
Konofagou EE, Provost J. Electromechanical wave imaging for noninvasive mapping of the 3D electrical activation sequence in caninesand humans in vivo. J Biomech 2012;45:856–864.
Marwick TH. Measurement of strain and strain rate by echocardiography: Ready for prime time?. J Am Coll Cardiol 2006;47:1313–
1327.
Mason SJ, Weiss JL, Weisfeldt ML, Garrison JB, Fortuin NJ. Exercise
echocardiography: Detection of wall motion abnormalities during
ischemia. Circulation 1979;59:50–59.
Matsuno Y, Taki H, Yamamoto H, Hirano M, Morosawa S, Shimokawa H, Kanai H. Ultrasound imaging of propagation of myocardial contraction for non-invasive identification of myocardialischemia. Jpn J Appl Phys 2017;56 07JF05.
McVeigh E, Faris O, Ennis D, Helm P, Evans F. Electromechanicalmapping with MRI tagging and epicardial sock electrodes. J Electro Cardiol 2002;35:61–64.
Pernot M, Fujikura K, Fung-Kee-Fung SD, Konofagou EE. ECGgated, mechanical and electromechanical wave imaging of cardiovascular tissues in vivo. Ultrasound Med Biol 2007;33:1075–1085.
Picano E, Lattanzi F, Orlandini A, Marini C, L’Abbate A. Stress echocardiography and the human factor: The importance of being
expert. J Am Coll Cardiol 1991;17:666–669.
Pislaru C, Bruce CJ, Anagnostopoulos PC, Allen JL, Seward JB, Pellikka PA, Ritman EL, Greenleaf JF. Ultrasound strain imaging ofaltered myocardial stiffness stunned versus infarcted reperfused
myocardium. Circulation 2004;109:2905–2910.
Ramanathan C, Ghanem RN, Jia P, Ryu K, Rudy Y. Noninvasive electrocardiographic imaging for cardiac electrophysiology andarrhythmia. Nat Med 2004;10:422–428.
Sabbah HN, Marzilli M, Stein PD. The relative role of subendocardium
and subepicardium in left ventricular mechanics. Am J Physiol
1981;240:H920–H926.

Volume 47, Number 4, 2021
Salles S, Chee AJY, Garcia D, Yu ACH, Vray D, Liebgott H. 2-D arterial wall motion imaging using ultrafast ultrasound and transverse
oscillations. IEEE Trans Ultrason Ferroelectr Freq Control
2015;62:1047–1058.
Sasayama S, Franklin D, Ross J, Jr, Kemper WS, McKown D.
Dynamic changes in left ventricular wall thickness and their use inanalyzing cardiac function in the conscious dog. Am J Cardiol
1976;38:870–879.
Shapiro E, Marier DL, St John Sutton MG, Gibson DG. Regional nonuniformity of wall dynamics in normal ventricle. Br Heart J
1981;45:264–270.
Sutherland GR, Salvo GD, Claus P, D’hooge J, Bijnens B. Strain andstrain rate imaging: A new clinical approach to quantifying
regional myocardial function. J Am Soc Echocardiogr 2004;17:
788–802.
Tanaka M, Sakamoto T, Sugawara S, Katahira Y, Tabuchi H, Nakajima H, Kurokawa T, Kanai H, Hasegawa H, Ohtsuki S. A newconcept of the contraction-extension property of the left ventricularmyocardium. J Cardiol 2014a;63:313–319.
Tanaka M, Sakamoto T, Sugawara S, Katahira Y, Tabuchi H, Nakajima H, Kurokawa T, Kanai H, Hasegawa H, Ohtsuki S. Non-uniform distribution of the contraction/extension (C-E) in the left
ventricular myocardium related to the myocardial function. J Cardiol 2014b;64:401–408.
Tanaka M, Sakamoto T, Sugawara S, Katahira Y, Hasegawa K,
Nakajima H, Kurokawa T, Kanai H, Hasegawa H. Deformability of the pulsating left ventricular wall: A new aspect elucidated by high resolution ultrasonic methods. J Cardiol
2017;69:462–470.
Tanaka M, Sakamoto T, Saijo Y, Katahira Y, Sugawara S, Nakajima H,
Kurokawa T, Kanai H. Role of intra‑ventricular vortex in left ventricular ejection elucidated by echo‑dynamography. J Med Ultrason
2019;46:413–423.
Theroux P, Franklin D, Ross J, Jr, Kemper WS. Regional myocardial
function during acute coronary artery occlusion and its modificationby pharmacologic agents in the dog. Circ Res 1974;35:896–908.
Yoshiara H, Hasegawa H, Kanai H. Ultrasonic imaging of propagationof contraction and relaxation in the heart walls at high temporal resolution. Jpn J Appl Phys 2007;46:4889–4896.