IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022

1353

Appropriate Window Function and Window
Length in Multifrequency Velocity Estimatorfor Rapid Motion and Localityof Layered Myocardium
Yu Obara , Shohei Mori , Member, IEEE, Mototaka Arakawa , Member, IEEE,and Hiroshi Kanai, Senior Member, IEEE

Abstract — The heart wall has a multilayered structure
and moves rapidly during ejection and rapid filling periods.
Local strain rate (SR) measurements of each myocardial
layer can contribute to accurate and sensitive evaluationsof myocardial function. However, ultrasound-based velocity estimators using a single-frequency phase difference
cannot realize these measurements owing to insufficientmaximum detectable velocity, which is limited by a quadrature frequency. We previously proposed a velocity estimator using multifrequency phase differences to improve themaximum detectable velocity. However, the improvement
is affected by a spatial discrete Fourier transform (DFT)window length that represents the locality of the velocity estimation. In this article, we theoretically describethat shortening the window increases the interference
between different frequency components and decreases
the maximum detectable velocity. The tradeoff between the
maximum detectable velocity and the window length was
confirmed through simulations and a water-tank experiment. Under the tradeoff, the Hanning window, which wasused in previous studies, is not always appropriate forthe local measurement of the velocity, which sometimes
exceeds 100 mm/s depending on the subject, direction of
the ultrasound beam to the heart wall, and cardiac periods.
In the in vivo measurement with the short window, the Tukey
window with a large flat part that has a high-frequency
resolution and ameliorates the discontinuity at both edges
Manuscript received December 27, 2021; accepted February 17, 2022.
Date of publication February 21, 2022; date of current version March 28,
2022. This work was supported in part by the Japan Society for the
Promotion of Science KAKENHI under Grant 20H02156 and Grant
19K22943, and in part by the Ministry of Education, Culture, Sports,
Science and Technology, University Fellowship Funding Project for
Innovation Creation in Science and Technology. (Corresponding author:
Yu Obara.)
This work involved human subjects or animals in its research. Approvalof all ethical and experimental procedures and protocols was granted by
the Ethics Committee of the Graduate School of Engineering, Tohoku
University, Sendai, Japan, under Approval No. 08A-2.
Yu Obara is with the Graduate School of Biomedical Engineering, Tohoku University, Sendai 9808579, Japan (e-mail:yu.obara.s2@dc.tohoku.ac.jp).
Shohei Mori is with the Graduate School of Engineering, Tohoku University, Sendai 9808579, Japan (e-mail: shohei.mori.a2@tohoku.ac.jp).
Mototaka Arakawa and Hiroshi Kanai are with the Graduate
School of Biomedical Engineering and the Graduate School ofEngineering, Tohoku University, Sendai 9808579, Japan (e-mail:mototaka.arakawa.c3@tohoku.ac.jp; hiroshi.kanai.e7@tohoku.ac.jp).
Digital Object Identifier 10.1109/TUFFC.2022.3153048
.

.

.

.

.

.

.

.

of the windowed signal was appropriate to measure themaximum velocity. This study offers the potential for local
measurements of each myocardial layer using the multifrequency velocity estimator with the appropriate windowfunction and window length.

Index Terms — Aliasing, discrete Fourier transform (DFT),
multifrequency phase differences, velocity measurement,
window function, window length.

I. I NTRODUCTIONM

YOCARDIAL ischemia leads to ischemic heart diseases, such as myocardial infarction and angina, which
are the main causes of death worldwide [1]. Myocardialischemia in the early stages can be repaired using reperfusion
therapy. It is, therefore, important to detect it at an early stage.
Because it occurs in the subendocardial region of the heart
wall and progresses transmurally toward the subepicardial
region [2]–[4], the local evaluation of myocardial function in
the transmural direction can be useful for its detection.
Ultrasound-based measurements of strain rate (SR) havebeen introduced to noninvasively evaluate regional myocardial
function [5]–[8]. The myocardial SR, in which the whole
motion of the heart wall is removed, represents an instantaneous change in the myocardial thickness. Previous studies
have suggested that myocardial SR can assess an acutely
ischemic myocardium [9]–[12]. Because the heart wall has
a multilayered structure, local SR measurements of eachmyocardial layer can contribute to the highly accurate andsensitive evaluations of myocardial function.
Moreover, the SR distribution has the potential to reveal
the mechanism of the transition process in the myocardium
from relaxation to contraction. In our previous studies,
the propagation of the contractile response caused by the
conduction of electrical excitation was measured using avelocity waveform [13]–[15] or an SR distribution of the
myocardium [16]–[18]. Because the whole motion of the heart
wall is removed in the SR measurement, it may be useful formeasuring the propagation of the local and minute contractile
responses [18], thereby contributing to the detailed elucidation

1525-8955 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See https://www.ieee.org/publications/rights/index.html for more information.
.

.

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

1354

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022of the transition process in the myocardium from relaxation
to contraction.
The myocardial SR distribution has been calculated from thevelocities at multiple positions in the heart wall. In addition,
the velocities in the heart wall are useful for understanding
whether the heart wall motion is related to the measured
SR. Thus, velocity estimation is important for SR measurement. A velocity-estimation method based on the phasedifferences between consecutive frames has been developed
for measuring the velocities in the heart wall [19], [20].
This method uses the phase differences between consecutive
frames in the demodulated signals and estimates the subpixel
displacements without interpolation. Thus, this method iscomputationally inexpensive compared to the block matching method [21] and has been applied for real-time SR
imaging [6], [22].
The maximum detectable velocity in this method, which isbased on the phase difference of a single frequency, is limited
and depends on the quadrature frequency. The maximum
detectable velocity may sometimes be insufficient to measure
the rapid motion of the heart wall, especially during theejection and rapid filling periods. In these cases, aliasing
occurs in the velocity waveform even when the frame rate
is set to several hundreds of Hz. Although the aliasing often
occurs in the normal myocardium rather than the ischemic
myocardium because the velocity decreases owing to the
myocardial ischemia, the aliasing is a problem to detect
the abnormality owing to the myocardial ischemia because
the large velocity in the normal myocardium has to be correctly measured as a reference.
To overcome this problem, several studies have proposed
anti-aliasing techniques based on the multifrequency phasedifferences between two consecutive frames [23]–[25]. The
improvement of the maximum detectable velocity in themeasurement of blood flow by a large spatial window has
already been confirmed [23]–[25]. To estimate the velocity
in the heart wall, we previously proposed a multifrequency
phased-tracking method [26] that achieves a larger maximum detectable velocity than that achievable using a singlefrequency phase difference.
The key factor that determines the maximum detectablevelocity in measurements based on the multifrequency phasedifferences is the length of the spatial window. The multifrequency phase differences should be measured independently
for each frequency component. However, a short spatial window that has a poor frequency resolution makes such measurements difficult. Thus, the improvement of the maximum
detectable velocity is negatively affected when the spatial
window for the velocity estimation is shortened.
By contrast, the local SR measurement in the heart wall
requires the local velocity to be estimated with a short spatial
window. The heart wall has a large velocity component related
to its whole motion and a minute velocity component causedby the local change in the thickness of the myocardial layer.
Because the velocity estimation assumes that the velocities are
uniform within the spatial window, it is difficult to measure theminute velocity caused by the local change in the thickness ofthe myocardial layer using a long spatial window. A previous

study confirmed that the spatial window length affects the SR
distribution in a layered myocardium [27].
There is a tradeoff between the maximum detectable velocity and the window length in a velocity estimation based on
the multifrequency phase differences. Thus, the influence of
the window length of the velocity estimation on the maximum
detectable velocity should be validated to determine the appropriate spatial window for the local SR measurement basedon the multifrequency phase differences. The purpose of thisstudy is to elucidate the relationship between the maximum
detectable velocity and the window length in the velocity
estimation. Therefore, we validate the relationship throughsimulations and a water-tank experiment. We then discuss
the tradeoff between the maximum detectable velocity andthe length of the spatial window and present an appropriate
window function for the rapid motion and locality of thelayered myocardium through an in vivo experiment.
II. P RINCIPLES

A. Velocity Estimation Based on Single-Frequency
Phase Difference
The displacement z(n) between the nth and the (n+n)thframes leads to a delay time between the received radio
frequency (RF) signals, s(t; n) and s(t; n + n), in the same
ultrasound beam. The delay time τ (n) is given by2 · z(n)(1)
c0where c0 is the speed of sound.
The demodulated signals are obtained by applying quadrature demodulation with the carrier frequency f to the RF
signals. The phase differenceθ f (n) between the demodulated
signals isτ (n) =

θ f (n) = −2π f · τ (n) = −4π f

z(n)
.
c0

(2)

In digital signal processing, this phase difference θ f (n)can be calculated from the cross-spectrum between the nth
and the (n + n)th frames [28]. The spectra S( f ; n) andS( f ; n + n) are obtained by applying discrete Fouriertransform (DFT) to the windowed signals. The cross-spectrumC( f ; n) between the frames is given byC( f ; n) = S ∗ ( f ; n) · S( f ; n + n)(3)where ∗ denotes complex conjugation.
The phase of the cross-spectrum  C( f ; n) is the phasedifference θ f (n) between the demodulated signals with thecarrier frequency f . In a previous study, it was confirmed thatthe phase of the cross-spectrum at the center frequency fc
should be used to estimate the minute velocity accurately [28].
Thus, the average velocity between frames is estimated using C( f c ; n) asv̂ fc (n) = −c0 f FR  C( fc ; n)
·
.
4π
fc

(4)

Here, f FR is the frame rate.
If the amplitude of the phase difference exceeds π, thatis, the delay time caused by the displacement exceeds half

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

OBARA et al.: APPROPRIATE WINDOW FUNCTION AND WINDOW LENGTH IN MULTIFREQUENCY VELOCITY ESTIMATOR

1355

Fig. 1. Schematic of the velocity estimation based on multifrequency phase differences. (a) Received RF signals and their demodulated signals at
the nth and the (n + Δn)th frames. (b) Frequency components of the spectra of the received RF signals. (c) Amplitude of the cross spectrum between
the nth and the (n + Δn)th frames. (d) Phase of the cross spectrum. Δf representing the frequency interval in the spectrum obtained by DFT.

the period of the center frequency f c , then the phase of
the cross-spectrum is wrapped and aliasing occurs. Aliasing
causes a large displacement error. By substituting π into C( f c ; n) in (4), the maximum detectable velocity |v̂ f ,MAX |c
z fc ,MAX | based on the single-frequency
and displacement |phase difference are given by

 

c0
v̂ f ,MAX  = z f c ,MAX  · f FR =· f FR .
(5)
c
4 fc

B. Velocity Estimation Based on Multifrequency Phase
Differences
Fig. 1 schematically illustrates the velocity estimation basedon multifrequency phase differences. Assuming that the displacement from the nth frame to the (n +n)th frame, z(n),is spatially uniform in the DFT window, the phase differenceis linearly proportional to the carrier frequency, according to(2). Thus, the gradient of the phase difference an is expressed
as
4π
(6)
an = − z(n).
c0This is estimated using the cross-spectrum phase. The phasegradient of the cross-spectrum is estimated using the weighted
least-squares method to ameliorate the negative effects of
the attenuation and interference of backscattered waves frommultiple scatterers in the myocardium in the velocity estimation [26]. The root mean squared error (RMSE) αn of thephase of the cross-spectrum, weighted by the amplitude of
the cross-spectrum, is defined with the arbitrary variable a asfollows:
  fs /2
 f =0 |C( f ; n)| · δ 2f (a)(7)αn (a) =  fs /2
f =0 |C( f ; n)|δ 2f (a) = min | C( f ; n) − a f + 2πm|2 .
(8)mHere, f s is the sampling frequency and m is an arbitrary
integer. Because the phase of the cross-spectrum is within therange of [−π, π], the phase a f must also fall within the samerange. When the value of a f falls outside this range, 2π is
subtracted from or added to the value of a f to handle thewrapped phase.
From (6), the average velocity v̂(n) between the frames isexpressed in terms of the estimated phase gradient ân as
c0 f FR· ân
4π
ân = arg min αn (a).

 · f FR = −v̂(n) = z
a

(9)
(10)

As shown in Fig. 1(d), the phase gradient can be estimatedbased on the multifrequency phase differences even when
the phase wraps at the center frequency. However, the phasegradient based on the single-frequency phase difference in(4),  C( f c ; n)/ f c , does not correspond to the true phasegradient. The multifrequency phase differences allow us tomeasure large velocities that cannot be measured using only
the single-frequency phase difference.
The phase gradient cannot be estimated when the change inthe phase of the cross-spectrum between adjacent frequency
components exceeds π. Thus, delay times that exceed half
the period of the DFT window in the time domain cannot
be detected from the cross-spectrum phase. The frequency
interval  f of the cross-spectrum phase is determined by thelength of the DFT window L DFT as follows:
f =

c0
.
2 · L DFT

(11)

Thus, the maximum phase gradient |â MAX |, displacement
 MAX |, and velocity |v̂ MAX | limited by the frequency interval|zAuthorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

1356

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022

 f are given bymodel
π
2π · L DFT
=
f
c0

(12)

L
z MAX | · f FR = DFT · f FR .
|v̂ MAX | = |2

(13)|â MAX | =This maximum velocity |v̂ MAX | is usually larger than the
maximum velocity |v̂ fc ,MAX | that can be detected based on
the single-frequency phase difference because the length of
the DFT window L DFT should be set to be longer than half
the wavelength at the center frequency f c [28].

C. Formulization of Multifrequency Phase Differences
Measured With DFT Window
In both the single-frequency and multifrequency phasedifference methods, it is assumed that the phase differenceat each frequency component is independent of the others,
that is, there is no interference between different frequency
components. When this assumption is correct, the maximum
detectable velocity based on the multifrequency phase differences corresponds to the maximum velocity |v̂ MAX | limited by
the frequency interval in (13).
In practice, however, there is interference between the
different frequency components because both the bandwidthof the received RF signal and the length of the window are
finite. Interference between the different frequency components should be considered, especially for local measurements.
Shortening the spatial DFT window improves the locality ofthe velocity estimation. However, it also results in poorer
frequency resolution and makes the independent measurementof each frequency component difficult. The assumption that
there is no interference between different frequency components then becomes incorrect. In this case, the maximum
detectable velocity based on the multifrequency phase differences becomes smaller than the maximum velocity |v̂ MAX |limited by the frequency interval in (13), as the DFT window
becomes shorter.
We formulized the measured cross-spectrum phase (i.e.,the multifrequency phase differences) obtained by DFT tocapture the influence of the interference between different
frequency components on the maximum detectable velocity
in the velocity estimation based on the multifrequency phasedifferences. In this section, a very simple model without
any noise is introduced to evaluate only the influence of
the interference between different frequency components on
the velocity estimation based on the multifrequency phasedifferences. The pulse wave representing the impulse responseh(t) of the ultrasound transducer is often approximated as the
product of a cosine wave at the center frequency f c with theGaussian function ash(t) = e−t /σ · cos(2π f c t).
2

2

(14)

Here, t = l/ f s (l = 0, 1, 2, . . .) is a discrete time and σ isthe standard deviation of the Gaussian function. When there
is a single point scatterer on the ultrasound beam, the received
RF signal sm (t; n) is expressed according to the convolutionsm (t; n) = A0 · h(t − tn )= A0 · e−(t−tn ) /σ · cos(2π f c (t − tn )).
2

2

(15)

Here, A0 is the reflectance of the scatterer and tn is thetwo-way propagation time from the ultrasound probe to thescatterer in the nth frame.
The true spectrum without the influence of windowing is
obtained by applying a discrete-time Fourier transform to theinfinite-length modeled RF signal. The true spectrum Sm ( f ; n)of the infinite-length modeled RF signal is given by the
convolution of the cosine wave at the center frequency f c and
the Gaussian function in the frequency domain as follows:
Sm ( f ; n) = A0 · e−σ ( f − fc ) · e− j 2π tn f
= |Sm ( f ; n)| · e− j 2π tn f .
2

2

(16)

From (3) and (16), the true cross-spectrum Cm ( f ; n) between
the true spectra Sm ( f ; n) and Sm ( f ; n + n) at the nth and(n + n)th frames, respectively, is expressed as
Cm ( f ; n) = A0 e−σ ( f − fc ) e j 2π tn f · A0 e−σ ( f − fc ) e− j 2π(tn +τ (n)) f= |Sm ( f ; n)|2 · e− j 2π τ (n) f .
(17)
2

2

2

2

Here, τ (n) is the delay time generated by the displacement from the nth frame to the (n + n)th frame. In thetrue cross-spectrum phase  Cm ( f ; n), only the delay timecomponent e− j 2π τ (n) f remains because the propagation timecomponent e− j 2π tn f cancels out. Thus,  Cm ( f ; n) is linearly
proportional to the frequency f , as seen in (2).
The spectrum obtained by DFT using a finite-length window
is expressed as the convolution of the true spectrum Sm ( f ; n)and the window function spectrum W ( f ). Thus, the measured
spectrum obtained by DFT, Sm ( f k ; n), is expressed asSm ( f k ; n) = (Sm ( f ; n) ∗ W ( f ))| f = fk
fk = k ·  f

(18)
(19)

where f k is a discrete frequency and k is an integer. The rangeof the discrete frequency f k is [− f s /2: fs /2]. In this article,
the integer k in (19) is defined as kc when f kc = f c .
Most of the window functions used for the DFT, such as
the Hanning window, are even functions. In this case, thephase  W0 ( f ) of the window function spectrum, in which thecenter point is set at the origin, has a zero-phase characteristic.
Thus, the spectrum of the window function W ( f ) has only
the linear phase component e− j 2π tn f , which is relevant to thecenter point of the window function set at the time tn , and thephase  W ( f ), which is proportional to the frequency f . The
spectrum of the window function W ( f ) is given by
W ( f ) = W0 ( f ) · e− j 2π tn f .

(20)

From (18)−(20), the measured spectrum at the nth frameSm ( f k ; n) that is obtained by DFT, is expressed as follows:Sm ( f k ; n) =+∞−∞Sm ( f ; n) · W ( f k − f )d f+∞

|Sm ( f ; n)| · W0 ( f k − f )d f
= e− j 2π tn fk−∞
 − j 2π t f
 
n k
.
(21)
=  Sm ( f k ; n) · eAuthorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

OBARA et al.: APPROPRIATE WINDOW FUNCTION AND WINDOW LENGTH IN MULTIFREQUENCY VELOCITY ESTIMATOR

1357

Fig. 2. (a) Amplitude of the spectrum of the modeled RF signal |Sm (f ; n)| and that of the window function spectrum W0 (f − fkc ) convoluted at thecenter frequency. (b) Phase of the true cross-spectrum  Cm (f ; n) and that of the measured cross-spectrum  Cm (fkc ; n) at the center frequency.
(c) Schematic of the summation of the delay time components e−j 2πτ (n)f in  Cm (fkc ; n). (d) Amplitude of the spectrum of the modeled RF signal|Sm (f; n)| and that of the window function spectrum W0 (f − fkc +1 ), convoluted at the adjacent frequency to the center frequency. (e) Phase of the true
cross-spectrum  Cm (f ; n) and that of the measured cross-spectrum  Cm (fkc +1 ; n) at the adjacent frequency to the center frequency. (f) Schematicof the summation of the delay time components e−j 2πτ (n)f in  Cm (fkc +1 ; n). Note the location of the mainlobe of the window function is not centeredwith the bandwidth of |Sm (f ; n)| when the frequency fk to be calculated does not correspond to the center frequency fc .

When the center point of the window function does not change
between the nth and the (n + n)th frames, the measured
spectrum at the (n + n)th frame Sm ( f k ; n + n), which is
obtained by DFT, is expressed as follows:Sm ( f k ; n + n) = e− j 2π tn f k

+∞
−∞

(22)

Cm ( fk ; n)=  Sm ( f k ; n) ·e+ j 2π tn fk=  S  ( f k ; n) ·m−∞|Sm ( f ; n)|·W0 ( fk − f )·e− j 2π τ (n) f d f+∞−∞

Cm ( fk ; n) = +∞−∞|Sm ( f ; n)|·W0 ( f k − f )·e− j 2π τ (n) f d f .
(24)

From (3) and (22), the measured cross-spectrum Cm ( f k ; n)between Sm ( f k ; n) and Sm ( f k ; n+n) is expressed as follows:+∞|Sm ( f ; n)|· W0 ( fk − f ) · e− j 2π τ (n) f d f.

· e− j 2π tn fk

The measured cross-spectrum phase  Cm ( f k ; n) is formulized
as follows:|Sm ( f ; n)|·W0 ( f k − f )·e− j 2π τ (n) f d f .
(23)Thus, the phase of the measured cross-spectrum  Cm ( f k ; n)corresponds to the phase of the sum of the delay timecomponent e− j 2π τ (n) f weighted by |Sm ( f ; n)| · W0 ( f k − f ).

D. Influence of the Window Length of the Velocity
Estimation Based on Multifrequency Phase Differenceson the Maximum Detectable Velocity
Fig. 2 schematically illustrates the influence of the interference between different frequency components on the measured
cross-spectrum phase used in the velocity estimation. Fig. 2(a)illustrates the amplitude of the true spectrum |Sm ( f ; n)| andthat of the window function W0 ( f − f kc ) shifted by the center
frequency. In Fig. 2(a), the mainlobe width corresponds tothe frequency resolution of the window function. Fig. 2(b)illustrates the phase of the true cross-spectrum  Cm ( f ; n) andthat of the measured cross-spectrum at the center frequency, C  ( f k ; n).
cmFig. 2(c) schematically illustrates the summation of theweighted delay time components |Sm ( f ; n)| · W0 ( f kc − f ) ·e− j 2π τ (n) f (blue vectors) at f k = f kc = fc , where both

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

1358

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022

Fig. 3 schematically illustrates the influence of the bias
errors in the measured cross-spectrum phases on the velocity estimation based on the multifrequency phase differences. Fig. 3(a) illustrates the amplitude of the cross-spectrum|Cm ( f ; n)| and that of the window functions W0 ( f − f kc ) andW0 ( f − f kc +1 ). Fig. 3(b) illustrates the wrapped true phasegradient, the unwrapped true phase gradient, the erroneously
estimated phase gradient, and the phase of the measured crossspectrum  Cm ( f k ; n).
When the measured cross-spectrum amplitude |Cm ( f k ; n)|is used as the weight to estimate the phase gradient, as shownin (7)–(10), the measured phase at the center frequency C  ( f k ; n) is dominant over those at other frequencies. Thecmmeasured cross-spectrum phase coincides with the true phaseat the center frequency. The possible values of the phasegradient an that passes through the origin and the phase at thecenter frequency are expressed asan (m) =Fig. 3. (a) Amplitude of the cross-spectrum |Cm (f ; n)| and that of
the window functions W0 (f − fkc ) and W0 (f − fkc +1 ). (b) Wrapped truephase gradient, the unwrapped true phase gradient, the erroneously estimated phase gradient, and the phase of the measured cross-spectrum, C (f ; n). Note that both the wrapped, and the unwrapped measured
m kcross spectrum  Cm (fk ; n) are illustrated.

|Sm ( f ; n)| and W0 ( f − f kc ) are symmetric with respect to thecenter frequency f c , as shown in black and blue in Fig. 2(a).
Thus, as shown in Fig. 2(c), the phase of the sum of theweighted delay time components,  (Cm ( f kc ; n)/|Sm ( f kc ; n)|),(white vector) corresponds to the true value  Cm ( fc ; n) at
the center frequency f c . The phase of the measured crossspectrum  Cm ( f kc ; n), therefore, coincides with the truephase  Cm ( f c ; n) in principle, as theoretically confirmed by
Hasegawa and Kanai [28].
Fig. 2(d) illustrates the amplitude of the true spectrum|Sm ( f ; n)| and that of the window function W0 ( f − f kc +1 )shifted by the frequency f kc +1 adjacent to the center frequency f kc . Fig. 2(e) illustrates the phase of the true crossspectrum  Cm ( f ; n) and that of the measured cross-spectrum C  ( f k +1 ; n). Fig. 2(f) schematically illustrates the summacmtion of the weighted delay time components |Sm ( f ; n)| ·W0 ( f kc +1 − f )·e− j 2π τ (n) f (red vectors), when f k = fkc +1 > f c .
W0 ( f − f kc +1 ) is symmetric with respect to the discrete
frequency f kc +1 , as shown in red in Fig. 2(d). However,|Sm ( f ; n)| is asymmetric with respect to the discrete frequency
f kc +1 , as shown in black in Fig. 2(d). Thus, a bias error occurs
in the phase of the measured cross-spectrum  Cm ( f kc +1 ; n),as shown in Fig. 2(e). As shown in Fig. 2(f), positive bias
errors occur when f k > f c . Conversely, the bias error values
are negative when f k < f c . The amount of bias error in thephase of the measured cross-spectrum depends mainly on themainlobe width of the window function W0 ( f − f k ) (i.e., thelength of the spatial window) compared with the bandwidth of
the true spectrum Sm ( f ; n) and the true phase gradient (i.e.,the delay time τ (n) generated by the displacement between
the frames). Because the velocity is estimated from the phasegradient, the amount of bias error also depends on the velocity.

Cm ( f c ; n) 2πm
+
.
fc
fc

(25)

Here, m is an arbitrary integer and the second term on the rightside of the equation accounts for the wrapped phase in the true
cross-spectrum phase at the center frequency,  Cm ( f c ; n). The
absolute value of m represents the number of phase wraps
in  Cm ( f k ; n) in the range of f k ≤ f c . For example, m =−1 means that the phase of the cross-spectrum is wrapped
once from a negative value to a positive value in the range off k ≤ fc , as shown in Fig. 3(b).
When the bias errors are large, the apparent gradient of
the measured phase of the cross-spectrum (red dashed line) issmaller than the true phase gradient (black lines) around thecenter frequency, as shown in Fig. 3(b). In this case, the bias
errors increase the weighted RMSE at the true phase gradientan (−1) and decrease that at an (0). Thus, the weighted RMSE
is minimum at an (0) (red solid line) even though the truephase gradient is an (−1), and an erroneous phase gradientis estimated. The large displacement error caused by thiserroneous estimation of m has the same effect as that causedby aliasing at the center frequency in the velocity estimation
based on a single-frequency phase difference. Therefore, themaximum detectable velocity becomes smaller than the maximum velocity limited by the frequency interval |v̂ MAX | in (13),owing to the bias errors caused by the interference between
different frequency components.

E. Influence of Noise on the Maximum Detectable
Velocity Based on Multifrequency Phase Differences
If noise presents, it causes a variance in the measured crossspectrum phase at each discrete frequency even if there is nointerference between different frequency components. Fig. 4schematically illustrates the influence of the noise variance inthe measured cross-spectrum phase on the velocity estimation
based on the multifrequency phase differences when the frequency resolution is sufficiently high. Fig. 4(a) illustrates theamplitude of the cross-spectrum |Cm ( f ; n)| and that of thewindow functions W0 ( f − f kc ) and W0 ( f − f kc +1 ). Fig. 4(b)illustrates the wrapped true phase gradient, the erroneously

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

OBARA et al.: APPROPRIATE WINDOW FUNCTION AND WINDOW LENGTH IN MULTIFREQUENCY VELOCITY ESTIMATOR

1359

frame rate. Therefore, the estimation error of the displacement
between frames was validated. The single-frequency phasedtracking method [20] or the multifrequency phased-tracking
method [26] was applied to the RF signals to estimate the displacement between the frames and the velocity. The detectablemaximum displacement was determined from the boundary atwhich m of the estimated phase gradient an (m) in (25) did
not coincide with the true value.
To present an appropriate window function for the rapid
motion and the locality of the layered myocardium, Tukey
windows were used. The Tukey window wT (t) is given by
Fig. 4. (a) Amplitude of the cross-spectrum |Cm (f ; n)| and that of thewindow functions W0 (f − fkc ) and W0 (f − fkc +1 ). (b) Wrapped true phasegradient, the wrapped erroneously estimated phase gradients an (0) andan (−2), and the phase of the measured cross-spectrum  Cm (fk ; n). Notethat there are only a few bias errors in the measured cross-spectrum
phases because the frequency resolution is sufficiently high.

estimated phase gradients, and the phase of the measured
cross-spectrum  Cm ( f k ; n).
If the noise variance is minimum at the center frequency
where the cross-spectrum amplitude |Cm ( f ; n)| is maximum,
its influence in the measured cross-spectrum phase on the possible phase gradients an (m) in (25) can be ameliorated by theweight of the measured cross-spectrum amplitude |Cm ( f k ; n)|.
However, the noise variance at a frequency other than thecenter frequency may increase the weighted RMSE at the truephase gradient an (−1) and decrease that at an (−2), as shownin Fig. 4(b). In this case, the noise variance in the measured
cross-spectrum phase causes the erroneous estimation of m inthe possible phase gradients an (m), and the phase gradient iserroneously estimated as an (−2). Thus, the large displacement
error of half the wavelength at the center frequency occurs
as an aliasing error. Depending on the noise variance, thephase gradient may be erroneously estimated as not onlyan (−2) but also an (0), an (1), an (−3), or the others withthe erroneous value of m. The large noise components may
cause an aliasing error even though the frequency resolutionof the DFT window is high enough to measure each frequency
component independently.
If the noise is regarded as white noise, the noise variancein the measured spectrum can be alleviated by elongating the
DFT window [29], [30]. In practice, however, the noise is
not only white noise in the in vivo measurement. Thus, theinfluence of the noise variance in the measured spectrum is
not always alleviated by elongating the DFT window.
III. M ATERIALS AND M ETHODSIn this study, the relationship between the maximum
detectable velocity and the window length was evaluatedby the estimation accuracy of m corresponding to the estimated phase gradient an (m) in (25) (i.e., the stability to thealiasing error). The erroneous estimation of m (i.e., aliasing error), owing to interference between different frequency
components, depends on the displacement between frames
which is determined by the velocity and the acquisitionwT (t) ⎧


1
2π 
r
⎪
⎪
1
+
cos
t
−
,
⎪
⎪
⎪
2
r
2
⎪
⎪
⎨
= 1,
⎪
⎪
⎪


⎪
⎪
r
1
2π 
⎪
⎪
t −1+
1 + cos
,
⎩
2
r
2

0≤t <

r
2

r
r
≤ t < 1−22r
1− ≤t < 12
(26)

where r is the parameter of the Tukey window. Dependingon the parameter r , the characteristic of the Tukey windowchanges between that of the Hanning window (r = 1) and therectangular window (r = 0).

A. Confirmation of the Principle
The modeled RF signal in (15) was used to confirm thebias errors in the measured cross-spectrum phases describedin Section II-D. The inverse of the standard deviation of theGaussian function in (15), σ −1 , was set to 1.5 MHz. The
delay time generated by the displacement between frames
was set to 0, 0.067, or 0.20 ms. The cross-spectrum wasmeasured with multiple window lengths. The length of thespatial window changes not only the frequency resolution but
also the frequency interval. Thus, zero points were added tothe end of the windowed signals in the time domain before
applying a DFT. This allowed the original frequency resolution
to be maintained [28].
B. Water-Tank Experiments
The influence of the window length of the velocity estimation on the maximum detectable velocity was validated in awater-tank experiment where the true displacement could bemeasured.
Fig. 5(a) shows a schematic of the water-tank experiment
system. A urethane rubber phantom was vibrated up and down
sinusoidally at 50 Hz, which was set based on the actual
velocity waveform of the heart wall [13], [20]. The maximum
velocity was set to approximately 40 mm/s.
The RF signals were acquired using an ultrasonic diagnostic apparatus (Prosound α-10; Hitachi-Aloka-Medical Ltd.,Tokyo, Japan) using a sector probe (UST-52101N; HitachiAloka-Medical Ltd., Tokyo, Japan) with a 3.75 MHz center
frequency. The transmitted wave was a plane wave withoutsteering and the received beam was formed without compounding. The sampling frequency was set to 15 MHz. As the

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

1360

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022

Fig. 5. (a) Schematic of the water tank experiment system. (b) Received RF signal at the first frame. (c) Velocity waveforms measured by laser
Doppler velocimetry or ultrasound measurement.

speed of sound in water was 1475 m/s and that in the phantom
was 1487 m/s, the RF inter-sample spacing in water was
0.049 mm and that in the phantom was 0.050 mm. The frame
rate of the ultrasound measurement was set to 1677 Hz.
The down-sampled frames were used as consecutive frames
to estimate the displacement of the phantom. Thus, the velocities of the phantom were estimated at multiple frame rates(1677, 839, 559, . . . , and 105 Hz) for the same data. The
maximum detectable velocity depends on the frame rate of
the ultrasound measurement. Thus, by decreasing the frame
rate in the velocity estimation, the displacement between theframes z(n) became larger than that at the original frame
rate. The time interval of the most down-sampled frames
at 105 Hz was similar to half the period of the sinusoidal
velocity waveform at 50 Hz. Because the phantom moved in
the same direction between the most down-sampled frames
at 105 Hz, the displacement between these frames had amaximum value of 0.25 mm.
A laser Doppler velocimeter (LV-1300; Ono Sokki, Kanagawa, Japan), which replaced the sector probe, was used forthe velocity measurement at the surface of the phantom. The
velocity measured by the laser Doppler velocimetry was usedas the reference. Fig. 5(c) illustrates the velocity waveforms
measured by the laser Doppler velocimetry or the ultrasoundmeasurement and the time window for evaluating the estimation error of the displacement between frames estimated by theultrasound measurement. As shown in Fig. 5(c), the trigger
signals measured by the ultrasound measurement without
down-sampling and laser Doppler velocimetry at different
timings were synchronized on a computer. The negative effectof the interference between different frequency components
depends on the delay time generated by the displacement
between frames. Therefore, the reference displacement z R (n)was obtained by integrating the velocity measured by the laser
Doppler velocimetry, v L (t), over a time duration corresponding

to the frame rate of the ultrasound measurement. The reference
displacement z R (n) is expressed as follows:z R (n) =(n+n)/ f FRv L (t)dt.

(27)

n/ f FRTo evaluate the estimation error, the time window was set to
ten times the wavelength of the velocity waveforms (200 ms),
as shown in Fig. 5(c). The estimation error of the displacement
between frames, e(n), was calculated for each pair of frames
as


(28)
e(n) = z(n) − z R (n).
The detectable maximum displacement was determinedfrom the boundary at which there was a large estimation error
e(n) caused by the erroneous estimation of m (i.e., aliasing
error).
The detectable maximum displacement was also obtained
from the simulation described in Section III-A. The inverse of
the standard deviation of the Gaussian function in (15), σ −1 ,was set to 3.0, 1.5, or 1.0 MHz. However, the envelope ofthe received RF signal in the water-tank experiment differedsignificantly from a Gaussian function and had the form of arectangular function, as shown in Fig. 5(b). Considering theactual waveform of the received RF signal, the received RF
signal at the first frame in Fig. 5(b) was used as the modeledRF signal at the nth frame sm (t; n) in (15).

C. Simulation for Determining the Appropriate Window
Function for the In Vivo Measurement
For determining the appropriate window function for the
rapid motion and the locality of the layered myocardium, thespeckle-like waveform in the heart wall and the noise in theactual in vivo measurement were simulated.
First, the speckle-like waveform in the heart wall was
simulated. Fig. 6(a) shows the schematic of the generation of

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

OBARA et al.: APPROPRIATE WINDOW FUNCTION AND WINDOW LENGTH IN MULTIFREQUENCY VELOCITY ESTIMATOR

1361

Fig. 6. (a) Schematics of the generation of the modeled RF signal in the simulation. (b) Modeled RF signals without noise, sm (t ; n) and sm (t ; n +Δn),in the simulation for the in vivo measurement. (c) Modeled RF signals with noise, sm,N (t ; n) and sm,N (t ; n + Δn), in the simulation for the in vivo
measurement.

the modeled RF signal in the simulation. In this study, themodeled RF signal was based on the convolution model of theimpulse response h(t) of the ultrasound transducer in (14) andthe scatterers’ distribution x(t; n). The inverse of the standard
deviation, σ −1 , was set to 2.1 MHz relating to the impulse
response of the ultrasound transducer that was obtained inmeasuring a 30 μm-diameter molybdenum wire with a plane
wave in a water tank. In simulating the in vivo measurement,
the multiple scatterers were distributed randomly in the beam
direction in space with a density of 15 scatterers per wavelength at the center frequency. The density was based on the
diameter of the myocardial fiber and its volume fraction in theheart wall [31]. The scatterers’ distribution at the (n + n)thframe was generated by adding an arbitrary delay time τ (n)to that at the nth frame, which represented the displacement
between the nth and the (n + n)th frames. Fig. 6(b) shows
the modeled RF signals, sm (t; n) and sm (t; n + n), in thesimulation of the in vivo measurement when the displacementfrom the nth frame to the (n+n)th frame was set at 0.14 mm.
As shown in Fig. 6(b), the envelopes of the modeled RF
signals, sm (t; n) and sm (t; n +n), were similar to the speckle
shown in in vivo measurements.
Next, the noise in the actual in vivo measurement was
simulated. In in vivo measurements, the change in the waveform of the RF signal occurs due to reasons other than justelectronic noise from an ultrasound diagnostic apparatus. Thisis because the movement in the orthogonal directions to thebeam direction and the nonuniform movement around theinterest point changes the scatterers’ distribution in the heartwall. The change in the waveform between the frames isregarded as a noise in the velocity estimation. Thus, this noise
was defined as the different components between the spectraS( f ; n) and S( f ; n+n) at the nth and the (n+n)th frames
other than the delay time component e− j 2π τ (n) f generated by
the displacement between the frames.
The noise spectrum N( f ; n) can be obtained by subtracting
the product of the spectrum at the nth frame S( f ; n) and thedelay time component e− j 2π τ (n) f from the spectrum at the(n + n)th frame S( f ; n + n). However, the true delay timeτ (n) cannot be known or measured in the in vivo measurement.
In this study, the estimated delay time τ̂ (n) using a long
DFT window (L DFT = 10 mm) was used to obtain the noise
spectrum. Because the long DFT window has a high-frequency
resolution, the spectra S( f ; n) and S( f ; n + n), and thedelay time component e− j 2π τ (n) f can be estimated without
the interference between the different frequency components.
The noise spectrum N( fk ; n) and the signal-to-noise ratio(SNR) ρ(n) in the velocity estimation were defined as
N( f k ; n) = S( f k ; n + n) − S( fk ; n) · e− j 2π τ̂ (n)· fk
 fs /2
2
f =0 S( f k ; n)ρ(n) = 10 log10  fsk/22
f k =0 N( f k ; n)

(29)
(30)

where S( fk ; n) and S( f k ; n + n) are the measured spectra at
the nth and the (n + n)th frames, respectively. Therefore,a sophisticated simulation can be realized using the noise
spectrum N( fk ; n) based on the in vivo measurement.
In this study, the noise spectrum was modeled based on thein vivo measurement to evaluate the average performance inthe velocity estimation in the heart wall. The noise spectrum

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

1362

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022

N( fk ; n) was obtained from the measurements of the interventricular septum (IVS) in the parasternal long-axis view
during a cardiac cycle, which will be explained in detail in
Section III-D. Because the bandwidth of the noise spectrum
was almost the same during a cardiac cycle, its amplitude
was generalized by the average over the cardiac cycle andwas used as that of the modeled noise spectrum |Nm ( f ; n)| inthe simulation. The phase of the noise spectrum  N( fk ; n)was different between adjacent frequency components andbetween the consecutive frames. Because the noise spectrum
phase was uniformly random in both the frequency domain
and the cardiac cycle, in this study, it was assumed that thephase of the noise spectrum was uniformly random. The phaseof the modeled noise spectrum  N m ( f ; n) was generalized
as a random value corresponding to the uniform distribution.
Considering the velocity-estimation SNR obtained in the invivo measurement, we added the modeled noise spectrum
Nm ( f ; n) to the spectrum of the modeled RF signal at the(n + n)th frame, Sm ( f ; n + n), as shown in Fig. 6(a). The
modeled RF signal at the nth frame may be affected by theelectronic noise from an ultrasound diagnostic apparatus withthe same power as that at the (n + n)th frame. Because theelectronic noise is often regarded as white noise, the white
noise component NW ( f ; n) included in the modeled noise
spectrum was added to the spectrum of the modeled RF signal
at the nth frame, Sm ( f ; n). The power of the white noise
component was almost a tenth of that of the modeled noise
spectrum. Fig. 6(c) shows the modeled RF signals with noise,sm,N (t; n) and sm,N (t; n + n), these are obtained by adding
the noise spectrum Nm ( f ; n) with the velocity-estimation SNRof 7 dB to the modeled RF signals without noise in Fig. 6(b).
As shown in Fig. 6(c), the envelope of the modeled RF signals
with noise at the (n + n)th frame, sm,N (t; n + n), wasdifferent from that at the nth frame, sm,N (t; n). This cannot be
expressed by adding only the white noise to the modeled RF
signal.
The relationship between the maximum detectable displacement and the window length was obtained using themodeled RF signals with noise, sm,N (t; n) and sm,N (t; n +n).
However, the detectable maximum displacement based on theestimation accuracy of m (i.e., the stability to the aliasing
error) is not uniquely determined because the scatterers’
distribution was random and the noise caused the variance inthe cross-spectrum phase. Thus, the probability of an aliasing
error was calculated for 100 generations of sm,N (t; n) andsm,N (t; n + n) for each displacement and window length.
The maximum detectable displacement was evaluated using
the calculated probability of the aliasing error.
In addition, the performance of the window function was
evaluated on the window length that can measure the maximum displacement between frames in the heart wall with thepermitted value of the probability of the aliasing error. In thisstudy, the permitted value of the probability of the aliasing
error was set to 0.05 or 0.10. Furthermore, in this study, themaximum velocity in the transmural direction of the heart wall
was obtained using the long DFT window (L DFT = 10 mm)
similar to the noise spectrum and the SNR in the velocity
estimation. The appropriate window function and window

length for the local velocity estimation was determined basedon the result of the simulation.

D. In Vivo Experiments
In the in vivo experiment, the improvement of the maximum detectable velocity by using the multifrequency phasedifferences was validated and the performance of the windowfunction in improving the maximum detectable velocity wascompared.
In vivo measurements were applied to the IVS of a healthy
23-year-old subject. The study was approved by the Ethics
Committee of Graduate School of Engineering, Tohoku University, Sendai, Japan, and the subject agreed to participate
in this study. The RF signal was acquired using an ultrasonic diagnostic apparatus (Prosound α-10; Hitachi-AlokaMedical Ltd., Tokyo, Japan) with a sector probe (UST52101N; Hitachi-Aloka-Medical Ltd., Tokyo, Japan) having
a 3.75 MHz center frequency. The seven plane waves inseven different directions at the steering angles of 6◦ intervals
were transmitted and compounded to obtain each received
beam [32]. The sampling frequency was set to 15 MHz.
Since the speed of sound in living tissue was assumed to
be 1540 m/s, the RF inter-sample spacing in the tissue wasestimated as 0.051 mm. The frame rate of the ultrasoundmeasurement was set at 860 Hz. The electrocardiogram (ECG)
was acquired using the three-point lead method, and the
ECG waveform was obtained from lead II. The phonocardiogram (PCG) waveform was measured using a small accelerometer attached to the subject’s chest.
The improvement of the maximum detectable velocity byusing the multifrequency phase differences was validated bycomparing the velocity waveform estimated locally by thesingle-frequency or multifrequency phased-tracking method.
The interest point of the velocity estimation was tracked in the
beam direction by integrating the estimated velocity waveform.
The performances of the window functions in improving the
maximum detectable velocity were compared by estimating
the velocities by applying different window functions with thesame length.
IV. R ESULTS

A. Bias Errors in the Measured Cross-Spectrum Phases
Fig. 7 shows the measured cross-spectrum between themodeled RF signals without noise obtained using multiple
window lengths when a single point scatterer was set on theultrasound beam. The color bar shows the length of the DFT
window. Fig. 7(a)–(c) shows the cross-spectrum at delay timesof 0, 0.067, and 0.20 ms, respectively.
In Fig. 7(a), where the delay time was set at 0 ms, themeasured phase corresponded to the true phase at all discrete
frequencies. This was because both the true cross-spectrum in(17) and the measured cross-spectrum in (23), in which 0 ms
was substituted for τ (n), had a zero-phase characteristic evenwhen the frequency resolution was poor (red points). Thus,
the influence of the interference between different frequency
components was small when the delay time was small.

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

OBARA et al.: APPROPRIATE WINDOW FUNCTION AND WINDOW LENGTH IN MULTIFREQUENCY VELOCITY ESTIMATOR

1363

Fig. 7. Measured cross-spectrum Cm (fk ; n) obtained using multiple window lengths. The cross-spectrum with delay times set at (a) 0, (b) 0.067,
and (c) 0.20 ms.

In Fig. 7(b), the true phase gradient was given by an (0) in(25) because the delay time (0.067 ms) did not exceed half
the period of the center frequency. The bias errors became
large and the apparent phase gradient decreased with thedecrease in frequency resolution (red points). However, thetrue phase gradient could be estimated because the true valueof m was minimal and the noise component was not added tothe modeled RF signal. Thus, the erroneous estimation of m inthe possible phase gradient did not occur when the phase of thecross-spectrum did not wrap and there was no additional noise.
In Fig. 7(c), the true phase gradient was given by an (−1)
in (25) because the phase of the cross-spectrum was wrapped
once from a negative value to a positive value. The bias errors
increased as the frequency resolution decreased (red points)
and as the apparent phase gradient decreased. These bias errors
caused the erroneous estimation of m in the possible phasegradient an (m) in (25). Thus, the phase gradient was estimatednot as an (−1) (purple solid line), but as an (0) (red solid line)
when the frequency resolution was poor (red points). Thiserroneous estimation resulted in a large displacement error.
Thus, in Fig. 7(c), the influence of the interference between
different frequency components on the velocity estimation
could not be suppressed despite the use of the cross-spectrum
amplitude as the weighting function.

B. Water-Tank ExperimentFig. 8 shows the spectra of the window function andthe estimation errors in the single-frequency and multifrequency phased-tracking methods in the water-tank experiment.
In Fig. 8(b)–(d), the horizontal axis represents the length of thewindow function and the vertical axis represents the reference
displacement between frames, z R . The color bar shows theestimation error of the displacement between frames, e(n).
The blue points represent a small estimation error, and the redpoints represent a large estimation error. The red dashed line
shows the maximum detectable displacement in the single-

 fc ,MAX |. The green
frequency phased-tracking method, |zz MAX | limited
solid line shows the maximum displacement |
by the frequency interval.
As shown in Fig. 8(b), a large displacement error occurredbecause of aliasing when the displacement of the phantom exceeded the maximum detectable displacement between
z fc ,MAX | < z R ). The boundary between the small
frames (|displacement errors (blue points) and the large displacement
errors (red points) is considered to be the maximum detectabledisplacement because the erroneous estimation of m occurred
there. The maximum detectable displacement in the singlefrequency phased-tracking method did not depend on thelength of the window function as is evident from (5).
As shown in Fig. 8(c) and (d), the estimation errors in
the multifrequency phased-tracking method depend on thelength of the window function. The maximum detectabledisplacement was increased by the use of multifrequency phasedifferences, but it decreased as the length of the window wasdecreased. This is consistent with the simulation results (blacksolid line). The maximum detectable displacement was smaller
than the maximum displacement limited by the frequency
 MAX | in both the simulations and the water-tank
interval |zexperiment. This is because the poor frequency resolution
made it difficult to measure large displacements.
The relationship between the maximum detectable displacement and the window length obtained in the water-tank
experiment did not correspond to those in the simulation
using the impulse response in (14), which had an envelope
similar to a Gaussian function. This was observed even though
it was confirmed that the maximum detectable displacement
decreased as the length of the window was decreased in boththe simulation and the water-tank experiment. In Fig. 8(c)and (d), the purple solid lines show the maximum detectabledisplacement obtained in the simulation using the received
RF signal from the water-tank experiment. This result corresponded well with the maximum detectable displacement

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

1364

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022

Fig. 8. (a) Spectrum of window function, (b) estimation error of the displacement between frames in the single-frequency phased-tracking method witha Hanning window, (c) that in the multifrequency phased-tracking method with a Hanning window, and (d) that in the multifrequency phased-tracking
method with a rectangular window.

obtained in the water-tank experiment, that is, the relationshipbetween the maximum detectable displacement and the locality of the velocity estimation also depended on the waveformof the received RF signal.
The relationships between the maximum detectable displacement and the length of the DFT window differed according to the window function. As shown in Fig. 8(d), a larger
maximum detectable displacement can be obtained using a
rectangular window than by using a Hanning window. Because
the mainlobe width of the rectangular window is half that of
the Hanning window, as shown in Fig. 8(a), the frequency
resolution of the former is twice that of the latter. Thus, themaximum detectable displacement of the rectangular window
approximately corresponded to that of the Hanning window
with twice the length of the spatial window.
As shown in Fig. 8(d), the maximum detectable displacement of the rectangular window did not monotonically increase
with the increase in the window length in both the simulations
and the water-tank experiment. This might be due to thespectral leakage caused by the discontinuity at both edges of
the rectangular window. The spectral leakage is caused by thediscontinuity at both edges of the window [33]. When thewindow length is not an integer multiple of the wavelength,
there is a discontinuity in the phase of the signal between

the edges of the window, and the spectral leakage increases.
In addition, the delay time of the signal causes a discontinuity
in the amplitude of the signal at both edges of the window
and increases the spectral leakage because the envelope amplitude of the modeled RF signal is not constant. The spectral
leakage caused by the discontinuity at both edges negatively
affected the maximum detectable displacement and resulted
in the non-monotonic relationship between the maximum
detectable displacement and the length of the rectangular
window. By contrast, the Hanning window is a window witha low sidelobe level [34]. Thus, the relationship between the
maximum detectable displacement and the length of the Hanning window depended mainly on the frequency resolution andits maximum detectable displacement monotonically increased
with the increase in the window length, as shown in Fig. 8(c).

C. Determining the Appropriate Window Function for the
In Vivo Measurement
In vivo measurement, the maximum velocity estimated using
the long DFT window of 10 mm was 116 mm/s. Thus,
we predicted that the maximum displacement that should be MAX |, was 0.14 mm. In simmeasured in the simulation, |zulating the in vivo measurement, the SNR in the velocity

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

OBARA et al.: APPROPRIATE WINDOW FUNCTION AND WINDOW LENGTH IN MULTIFREQUENCY VELOCITY ESTIMATOR

1365

Fig. 9. (a) Spectrum of the Tukey window, probability of an aliasing error in the multifrequency phased-tracking method with the Tukey window
where the parameter (b) r = 1, (c) r = 0.50, and (d) r = 0.

estimation was set to 7 dB, which was obtained as in (30)using the long DFT window at the time when the velocity inthe heart wall was larger than the maximum detectable velocity
based on a single-frequency phase difference.
Fig. 9 shows the spectra of the window functions andthe probability of an aliasing error in the multifrequency
phased-tracking method in the simulation of the in vivo measurement. In Fig. 9(b)–(d) the horizontal axis represents thelength of the window function and the vertical axis represents
the true displacement between frames, z. The color bar
shows the probability of an aliasing error in the multifrequency
phased-tracking method. The blue area represents a stableestimation. The red area represents an unstable estimation
owing to the aliasing error. The red dashed line shows the
maximum detectable displacement in the single-frequency
 fc ,MAX |.
phased-tracking method, |z
As shown in Fig. 9(b) and (d), the area where the probabilityof an aliasing error was high (red area) in the result using
the Hanning window was broader than that in the result
using the rectangular window. This might be representativeof the difference in the frequency resolution depending on
the mainlobe width, which corresponded to the result of
the water-tank experiment in Fig. 8. However, as shown in
Fig. 9(d), it was not easy to decrease the probability of analiasing error when using the rectangular window even though
the window length was increased. This may be a negative
effect of the discontinuity at both edges of the windowed
signal. This negative effect was not ameliorated by elongating

the DFT window and could not be ignored in simulating thein vivo measurement.
As shown in Fig. 9(c), this negative effect was ameliorated
in the result using the Tukey window (r = 0.50). Because the
Tukey window (r = 0.50) has a higher frequency resolution
than the Hanning window as shown in Fig. 9(a), the maximum
detectable displacement of the Tukey window (r = 0.50) may
be larger than that of the Hanning window. There is a tradeoff
between the frequency resolution and the negative effect of the
discontinuity at both edges of the windowed signal. Dependingon the parameter r in (26), the characteristic of the Tukey
window changes between that of the Hanning window and therectangular window, as shown in Fig. 9(a). Thus, the Tukey
window is useful to optimize this tradeoff. The performance of
the Tukey window was evaluated for each value of parameterr in (26).
Fig. 10 shows the minimum value of the window lengththat can measure the maximum displacement between frames
in the heart wall (0.14 mm) with the permitted value of theprobability of an aliasing error. In Fig. 10, the permitted valueof the probability of the aliasing error was set to 0.05 (blackline) or 0.10 (blue line). As shown in Fig. 10, the Tukey
window (r = 0.30) shows the best performance to locallymeasure the maximum displacement of 0.14 mm in both thepermitted values of the probability of the aliasing error. Whenthe permitted value of the probability of the aliasing error was
set to 0.05, the shortest Tukey window (r = 0.30) was 2.5 mm.
In this window function, the length of the tapered part toward

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

1366

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022

Fig. 10. Minimum value of the Tukey window length that can measure
the maximum displacement between frames in the heart wall with thepermitted value of the probability of an aliasing error.

the edge approximately corresponded to two wavelengths(0.40 mm) at each edge. This might ameliorate the negative
effect of the discontinuity at both edges of the windowed signal
in the rectangular window.
Interestingly, the Hanning window, which has a low sidelobe
level, was inferior to the Tukey window (r = 0.30) evenwhen the noise was present. As shown in Fig. 6(a), the powerof the white noise component was almost a tenth of that of
the modeled noise spectrum Nm ( f ; n). Thus, the advantageof the low sidelobe level of the Hanning window might not
be effective. The Tukey window (r = 0.30) has a highfrequency resolution and ameliorates the negative effect ofthe discontinuity at both edges because the Tukey window(r = 0.30) has a flat part around the center and the tapered
part toward the edges, as in (26). Thus, the Tukey window(r = 0.30) shows a better performance in improving the
maximum detectable displacement than the Hanning window.

the velocity waveform estimated using the single-frequency
phased-tracking method with the 2.5 mm Hanning window(red line) or the 2.5 mm Tukey window (r = 0.30) (blueline). In the single-frequency phased-tracking method, there
was no significant difference in the number of aliasing errors
between the velocity waveforms estimated with the Hanningwindow and the Tukey window (r = 0.30). This was because
the maximum detectable displacement in the single-frequency
phased-tracking method does not depend on the frequency
resolution, as confirmed in the water-tank experiment.
As shown in Fig. 11(e), the number of aliasing errors in
the velocity waveforms decreased using the multifrequency
phased-tracking method. Thus, the maximum detectable velocity was improved by the use of the multifrequency phase differences. However, in the result of the Hanning window [red line
in Fig. 11(e)], aliasing errors were not completely suppressed,
and they also occurred around the time from the ECG R-wave
at 400 ms even though the velocity waveform was estimatedusing the multifrequency phase differences. In this case, thelength of the Hanning window should be increased to suppress
the aliasing errors. By contrast, in the result of the Tukey
window (r = 0.30) [blue line in Fig. 11(e)], an aliasing error
did not occur during a cardiac cycle due to the advantage ofa high-frequency resolution. In addition, at the time from the
ECG R-wave of 30 ms, when the estimated velocity using thelong DFT window of 10 mm was maximum, the velocities at
all the spatial points in the IVS (13674 points) were estimatedusing the multifrequency phased-tracking method with the
Hanning window or the Tukey window (r = 0.30). The
number of points where an aliasing error occurred is 608(4%) or 302 (2%) in using the Hanning window or the Tukey
window (r = 0.30), respectively. These results suggested that
the Tukey window (r = 0.30) was useful for improving the
maximum detectable velocity in the in vivo measurement while
maintaining the locality of the velocity estimation.
V. D ISCUSSION

D. In Vivo Experiment Using the Appropriate Window
Function Determined in the Simulation
In the result of the simulation considering the in vivo
measurement condition, the appropriate window function was
the Tukey window (r = 0.30), as shown in Fig. 10. In thein vivo experiment, the performance of the Tukey window(r = 0.30) in improving the maximum detectable velocity
was compared with that of the Hanning window, which
was used in our previous studies [18], [26], [27]. Fig. 11
shows the B-mode image of the IVS from the parasternal
long-axis view, the ECG waveform, the PCG waveform, andthe velocity waveforms estimated using the single-frequency
and multifrequency phased-tracking methods. In Fig. 11(a),the yellow dot represents the interest point of the velocity
estimation where the center of the window function wasset. In Fig. 11(d) and (e), the red dashed line represents thedetectable maximum velocity based on the single-frequency
phase difference, |v̂ fc ,MAX |.
As shown in Fig. 11(d), around the time from the ECG
R-wave of 30 and 640 ms, aliasing errors occurred in

The maximum detectable velocity was improved by using
the multifrequency phase difference in both the water-tank
experiment and the in vivo experiment. Even though this maximum detectable velocity was negatively affected by the interference between the different frequency components, it waslarger than the maximum detectable velocity based on a singlefrequency phase difference when the cross-spectrum amplitude
was used as the weighting function. Thus, the multifrequency
velocity estimator was useful for measuring the rapid motionof the heart wall.
In this study, the plane wave was transmitted from thesector probe. Diverging waves were used to measure thelarge field of view, and the measurement using diverging
waves was proposed for the strain measurement of the heart
wall [32], [35], [36]. There is no significant difference in thespatial resolution in the beam direction between the plane
wave and the diverging wave [32], [37]. Thus, the condition of
the transmitted wave may not affect the relationship between
the maximum detectable displacement and the window length.
However, the SNR in the velocity estimation may depend

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

OBARA et al.: APPROPRIATE WINDOW FUNCTION AND WINDOW LENGTH IN MULTIFREQUENCY VELOCITY ESTIMATOR

1367

Fig. 11. (a) B-mode image of the IVS of a healthy young subject in measurement from the parasternal long-axis view. (b) ECG waveform, (c) PCG
waveform, and (d) velocity waveform estimated using a single-frequency phased-tracking method with a 2.5 mm Hanning window or a 2.5 mm Tukey
window (r = 0.30). (e) Velocity waveform estimated using a multifrequency phased-tracking method with a 2.5 mm Hanning window or a 2.5 mm
Tukey window (r = 0.30). (LV: left ventricular, RV: right ventricular).

on the transmitted wave condition owing to the differencein the temporal resolution and the spatial resolution in theorthogonal direction to the beam direction. Because the noisein the velocity estimation is related to the movement of theheart wall, as described in Section III-C, the changes in thetemporal resolution and the spatial resolution in the orthogonal
directions may affect the power of the noise spectrum. Even
though the transmitted wave was not a plane wave, the SNR in
the velocity estimation can be validated by the present method.
Therefore, the methods of this study may be beneficial for the
SR measurement using diverging waves.
Let us discuss the determination of the appropriate window length based on the relationship between the maximum
detectable velocity and the window length for the local SR
measurement in the heart wall. In this study, the in vivo
experiment was applied to the IVS from the parasternal
long-axis view and the velocity in the transmural direction
was measured. The multifrequency phased-tracking method
estimated only the velocity in the beam direction. Thus, theheart wall was measured from the parasternal long-axis view,
short-axis view, or the apical view depending on the velocity
component that we want to measure. The maximum velocity
in the heart wall can be approximately predicted based on
the reference value for each measurement view, region of
the heart wall, and cardiac periods. The reference value ismeasurable by the multifrequency phased-tracking methodwith a long DFT window, which can detect a sufficiently highdisplacement, as described in Section III-C. The maximum
velocity in the longitudinal direction is larger than that in thetransmural direction [38], [39]; the appropriate window lengthfor estimating the longitudinal component may be longer thanthat for estimating the transmural one.
The 1-D strain and SR imaging method have been used
for clinical applications [7], [40]–[42], and the 2-D or 3-D
strain and SR imaging have also been studied [43]–[45]. The
multifrequency phased-tracking method can be expanded to

the 2-D or 3-D velocity estimator if needed [46], [47]. In the
beam direction of the 2-D or 3-D velocity estimator based on
the multifrequency phase differences, the relationship between
the maximum detectable velocity and the window length may
correspond to that in this study.
The computational cost in using the multifrequency phasedifferences increases as the number of the analysis discrete
frequencies used for estimating the cross-spectrum phase gradient increases. Thus, the computational cost of the velocity
estimator based on the multifrequency phase differences ishigher than that based on a single-frequency phase difference.
One of the methods to decrease the computational cost is to
impose a limitation on the frequency components used for
estimating the cross-spectrum phase gradient. It is our future
work to implement the multifrequency phased-tracking method
to an ultrasound diagnostic apparatus for real-time imaging.
VI. C ONCLUSION
In this article, the relationship between the maximum
detectable velocity and the window length in the velocity
estimation based on the multifrequency phase differences was
described theoretically. The interference between different frequency components in a short DFT window causes bias errors
in the phase differences at discrete frequencies other thanthe center frequency. This decreased the maximum detectablevelocity and caused a large estimation error in the velocity.
The tradeoff between the maximum detectable velocity andthe window length was confirmed through simulations anda water-tank experiment. Under the tradeoff, the Hanningwindow, which was used in previous studies, is not always
appropriate for the local measurement of the velocity, whichsometimes exceeds 100 mm/s depending on the subject, direction of the ultrasound beam to the heart wall, and cardiac
periods. In the in vivo measurement with the short window,
the Tukey window with a large flat part, which has a highfrequency resolution and ameliorates the discontinuity at both

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

1368

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 69, NO. 4, APRIL 2022

edges of the windowed signal, was appropriate for measuring
the maximum velocity of the heart wall. This study offers
the potential for the local velocity and SR measurementsof each myocardial layer using the multifrequency velocity
estimator with the appropriate window function and window
length.
R EFERENCES[1] K. J. Foreman et al., “Forecasting life expectancy, years of life lost,and all-cause and cause-specific mortality for 250 causes of death:
Reference and alternative scenarios for 2016–40 for 195 countries andterritories,” Lancet, vol. 392, no. 10159, pp. 2052–2090, Nov. 2018, doi:
10.1016/S0140-6736(18)31694-5.
[2] K. A. Reimer, J. E. Lowe, M. M. Rasmussen, and R. B. Jennings, “Thewavefront phenomenon of ischemic cell death. 1. myocardial infarct sizevs duration of coronary occlusion in dogs,” Circulation, vol. 56, no. 5,
pp. 786–794, Nov. 1977, doi: 10.1161/01.CIR.56.5.786.
[3] H. Fujiwara, M. Ashraf, S. Sato, and R. W. Millard, “Transmural
cellular damage and blood flow distribution in early ischemia in pig
hearts,” Circulat. Res., vol. 51, no. 6, pp. 683–693, Dec. 1982, doi:
10.1161/01.RES.51.6.683.
[4] H. Kanaide, Y. Taira, and M. Nakamura, “Transmural anoxic wave front
and regional dysfunction during early ischemia,” Amer. J. Physiol.-Heart
Circulatory Physiol., vol. 253, no. 2, pp. H240–H247, Aug. 1987, doi:
10.1152/ajpheart.1987.253.2.H240.
[5] H. Kanai, H. Hasegawa, N. Chubachi, Y. Koiwa, and M. Tanaka,“Noninvasive evaluation of local myocardial thickening and its colorcoded imaging,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control,
vol. 44, no. 4, pp. 752–768, Jul. 1997, doi: 10.1109/58.655190.
[6] A. Heimdal, A. Støylen, H. Torp, and T. Skjærpe, “Real-time strain rate
imaging of the left ventricle by ultrasound,” J. Amer. Soc. Echocardiography, vol. 11, no. 11, pp. 1013–1019, Nov. 1998, doi: 10.1016/S08947317(98)70151-8.
[7] J. D’hooge et al., “Regional strain and strain rate measurements bycardiac ultrasound: Principles, implementation and limitations,” Eur.
Heart J.-Cardiovascular Imag., vol. 1, no. 3, pp. 154–170, Sep. 2000,
doi: 10.1053/euje.2000.0031.
[8] R. B. Shenouda, I. Bytyçi, M. Sobhy, and M. Y. Henein, “Reduced
regional strain rate is the most accurate dysfunction in predicting
culprit lesions in patients with acute coronary syndrome,” Clin. Physiol.
Funct. Imag., vol. 40, no. 1, pp. 21–29, Jan. 2020, doi: 10.1111/cpf.
12597.
[9] C. Pislaru, M. Belohlavek, R. Y. Bae, T. P. Abraham, J. F. Greenleaf, andJ. B. Seward, “Regional asynchrony during acute myocardial ischemia
quantified by ultrasound strain rate imaging,” J. Amer. College Cardiol.,vol. 37, no. 4, pp. 1141–1148, Mar. 2001.
[10] M. Belohlavek, C. Pislaru, R. Y. Bae, J. F. Greenleaf, and J. B. Seward,“Real-time strain rate echocardiographic imaging: Temporal and spatial
analysis of postsystolic compression in acutely ischemic myocardium,”
J. Amer. Soc. Echocardiography, vol. 14, no. 5, pp. 360–369, May 2001,
doi: 10.1067/mje.2001.110786.
[11] C. Pislaru, P. C. Anagnostopoulos, J. B. Seward, J. F. Greenleaf,
and M. Belohlavek, “Higher myocardial strain rates during isovolumic
relaxation phase than during ejection characterize acutely ischemic
myocardium,” J. Amer. College Cardiol., vol. 40, no. 8, pp. 1487–1494,
Oct. 2002.
[12] F. Jamal et al., “Can changes in systolic longitudinal deformation
quantify regional myocardial function after an acute infarction? An
ultrasonic strain rate and strain study,” J. Amer. Soc. Echocardiography, vol. 15, no. 7, pp. 723–730, Jul. 2002, doi: 10.1067/mje.2002.
118913.
[13] H. Kanai, “Propagation of vibration caused by electrical excitationin the normal human heart,” Ultrasound Med. Biol., vol. 35, no. 6,
pp. 936–948, Jun. 2009, doi: 10.1016/j.ultrasmedbio.2008.12.013.
[14] A. Hayashi, S. Mori, M. Arakawa, and H. Kanai, “Local twodimensional distribution of propagation speed of myocardial contraction
for ultrasonic visualization of contraction propagation,” Jpn. J. Appl.
Phys., vol. 58, no. SG, Jun. 2019, Art. no. SGGE05, doi: 10.7567/13474065/ab0d0b.
[15] H. Kanai, S. Mori, and M. Arakawa, “Speeds of contraction
responses propagating along septum at pre-ejection period are different between radial and longitudinal directions,” in Proc. IEEE
Int. Ultrason. Symp. (IUS), Oct. 2019, pp. 1418–1420, doi: 10.1109/
ULTSYM.2019.8925990.

[16] H. Yoshiara, H. Hasegawa, H. Kanai, and M. Tanaka, “Ultrasonicimaging of propagation of contraction and relaxation in the heart walls
at high temporal resolution,” Jpn. J. Appl. Phys., vol. 46, no. 7B,
pp. 4889–4896, Jul. 2007, doi: 10.1143/JJAP.46.4889.
[17] M. Tanaka et al., “A new concept of the contraction–extension property of the left ventricular myocardium,” J. Cardiol., vol. 63, no. 4,
pp. 313–319, Apr. 2014, doi: 10.1016/j.jjcc.2013.09.009.
[18] Y. Obara, S. Mori, M. Arakawa, and H. Kanai, “Measurement of propagation of local and minute contractile response in layered myocardium,”
Jpn. J. Appl. Phys., vol. 60, Mar. 2021, Art. no. SDDE02, doi: 10.35848/
1347-4065/abeabf.
[19] W. N. McDicken, G. R. Sutherland, C. M. Moran, and L. N. Gordon,“Color Doppler velocity imaging of the myocardium,” UltrasoundMed. Biol., vol. 18, nos. 6–7, pp. 651–654, 1992, doi: 10.1016/03015629(92)90080-T.
[20] H. Kanai, M. Sato, Y. Koiwa, and N. Chubachi, “Transcutaneous
measurement and spectrum analysis of heart wall vibrations,” IEEE
Trans. Ultrason., Ferroelectr., Freq. Control, vol. 43, no. 5, pp. 791–810,
Sep. 1996, doi: 10.1109/58.535480.
[21] C. Jia et al., “Comparison of 2-D speckle tracking and tissue Doppler
imaging in an isolated rabbit heart model,” IEEE Trans. Ultrason.,Ferroelectr., Freq. Control, vol. 57, no. 11, pp. 2491–2502, Nov. 2010,
doi: 10.1109/TUFFC.2010.1715.
[22] H. Kanai, Y. Koiwa, and J. Zhang, “Real-time measurements oflocal myocardium motion and arterial wall thickening,” IEEE Trans.
Ultrason., Ferroelectr., Freq. Control, vol. 46, no. 5, pp. 1229–1241,
Sep. 1999, doi: 10.1109/58.796128.
[23] H. J. Nitzpon, J. C. Rajaonah, C. B. Burckhardt, B. Dousse, andJ.-J. Meister, “A new pulsed wave Doppler ultrasound system to measure
blood velocities beyond the Nyquist limit,” IEEE Trans. Ultrason.,Ferroelectr., Freq. Control, vol. 42, no. 2, pp. 265–279, Mar. 1995, doi:
10.1109/58.365240.
[24] M. Mozumi, R. Nagaoka, M. Cinthio, and H. Hasegawa, “Anti-aliasing
method for ultrasonic 2D phase-sensitive motion estimator,” Jpn. J. Appl.
Phys., vol. 59, no. SK, May 2020, Art. no. SKKE20, doi: 10.35848/13474065/ab8c9a.
[25] J. Poree et al., “Dealiasing high-frame-rate color Doppler using dualwavelength processing,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control, vol. 68, no. 6, pp. 2117–2128, Jun. 2021, doi: 10.1109/TUFFC.
2021.3056932.
[26] Y. Obara, S. Mori, M. Arakawa, and H. Kanai, “Multifrequency phased
tracking method for estimating velocity in heart wall,” UltrasoundMed. Biol., vol. 47, no. 4, pp. 1077–1088, Apr. 2021, doi: 10.1016/j.
ultrasmedbio.2020.12.011.
[27] Y. Obara, S. Mori, M. Arakawa, and H. Kanai, “Strain rate distribution
in layered myocardium measured using local velocity estimator withmultifrequency phase differences,” Ultrasound Med. Biol., vol. 47, no. 9,
pp. 2768–2773, Sep. 2021, doi: 10.1016/j.ultrasmedbio.2021.05.021.
[28] H. Hasegawa and H. Kanai, “Modification of the phased-tracking
method for reduction of artifacts in estimated artery wall deformation,”
IEEE Trans. Ultrason., Ferroelectr., Freq. Control, vol. 53, no. 11,
pp. 2050–2064, Nov. 2006, doi: 10.1109/TUFFC.2006.145.
[29] A. P. G. Hoeks, T. G. J. Arts, P. J. Brands, and R. S. Reneman,“Comparison of the performance of the RF cross correlation and Doppler
autocorrelation technique to estimate the mean velocity of simulated
ultrasound signals,” Ultrasound Med. Biol., vol. 19, no. 9, pp. 727–740,
1993, doi: 10.1016/0301-5629(93)90090-B.
[30] Y. Honjo, H. Hasegawa, and H. Kanai, “Optimization of correlation
kernel size for accurate estimation of myocardial contraction and relaxation,” Jpn. J. Appl. Phys., vol. 51, Jul. 2012, Art. no. 07GF06, doi:
10.1143/JJAP.51.07GF06.
[31] I. J. LeGrice, B. H. Smaill, L. Z. Chai, S. G. Edgar, J. B. Gavin,and P. J. Hunter, “Laminar structure of the heart: Ventricular myocyte
arrangement and connective tissue architecture in the dog,” Amer. J.
Physiol.-Heart Circulatory Physiol., vol. 269, no. 2, pp. H571–H582,
Aug. 1995, doi: 10.1152/ajpheart.1995.269.2.H571.
[32] H. Hasegawa and H. Kanai, “High-frame-rate echocardiography using
diverging transmit beams and parallel receive beamforming,” J. Med.
Ultrason., vol. 38, no. 3, pp. 129–140, 2011, doi: 10.1007/s10396-0110304-0.
[33] F. J. Harris, “On the use of windows for harmonic analysis with thediscrete Fourier transform,” Proc. IEEE, vol. 66, no. 1, pp. 51–83,
Jan. 1978, doi: 10.1109/PROC.1978.10837.
[34] A. H. Nuttall, “Some windows with very good sidelobe behavior,” IEEE
Trans. Acoust., Speech, Signal Process., vol. ASSP-29, no. 1, pp. 84–91,
Feb. 1981, doi: 10.1109/TASSP.1981.1163506.

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.

OBARA et al.: APPROPRIATE WINDOW FUNCTION AND WINDOW LENGTH IN MULTIFREQUENCY VELOCITY ESTIMATOR[35] C. Papadacci, M. Pernot, M. Couade, M. Fink, and M. Tanter, “Highcontrast ultrafast imaging of the heart,” IEEE Trans. Ultrason., Ferroelectr., Freq. Control, vol. 61, no. 2, pp. 288–301, Feb. 2014, doi:
10.1109/TUFFC.2014.6722614.
[36] J. Grondin, V. Sayseng, and E. E. Konofagou, “Cardiac strain imaging
with coherent compounding of diverging waves,” IEEE Trans. Ultrason.,Ferroelectr., Freq. Control, vol. 64, no. 8, pp. 1212–1222, Aug. 2017,
doi: 10.1109/TUFFC.2017.2717792.
[37] L. Tong, H. Gao, H. F. Choi, and J. D’hooge, “Comparison of conventional parallel beamforming with plane wave and diverging wave
imaging for cardiac applications: A simulation study,” IEEE Trans.
Ultrason., Ferroelectr., Freq. Control, vol. 59, no. 8, pp. 1654–1663,
Aug. 2012, doi: 10.1109/TUFFC.2012.2371.
[38] P. Palka, A. Lange, A. D. Fleming, G. R. Sutherland, L. N. Fenn, andW. N. McDicken, “Doppler tissue imaging: Myocardial wall motion
velocities in normal subjects,” J. Amer. Soc. Echocardiogr., vol. 8, no. 5,
pp. 659–668, Sep.-Oct. 1995, doi: 10.1016/S0894-7317(05)80380-3.
[39] H. Dalen, A. Thorstensen, L. J. Vatten, S. A. Aase, and A. Stoylen,“Reference values and distribution of conventional echocardiographic
Doppler measures and longitudinal tissue Doppler velocities in a
population free from cardiovascular disease,” Circulation, Cardiovascular Imag., vol. 3, no. 5, pp. 614–622, Sep. 2010, doi: 10.1161/
CIRCIMAGING.109.926022.
[40] J. Koyama, P. A. Ray-Sequin, and R. H. Falk, “Longitudinal myocardial
function assessed by tissue velocity, strain, and strain rate tissue Doppler
echocardiography in patients with AL (primary) cardiac amyloidosis,” Circulation, vol. 107, no. 19, pp. 2446–2452, May 2003, doi:
10.1161/01.CIR.0000068313.67758.4F.
[41] J. Grondin, D. Wang, C. S. Grubb, N. Trayanova, and E. E. Konofagou,“4D cardiac electromechanical activation imaging,” Comput. Biol. Med.,
vol. 113, Oct. 2019, Art. no. 103382, doi: 10.1016/j.compbiomed.
2019.103382.
[42] V. Sayseng et al., “Catheter ablation lesion visualization with intracardiac strain imaging in canines and humans,” IEEE Trans. Ultrason.,Ferroelectr., Freq. Control, vol. 67, no. 9, pp. 1800–1810, Sep. 2020,
doi: 10.1109/TUFFC.2020.2987480.
[43] Y. Notomi et al., “Measurement of ventricular torsion by twodimensional ultrasound speckle tracking imaging,” J. Amer. College
Cardiol., vol. 45, no. 12, pp. 2034–2041, Jun. 2005.
[44] S. Langeland et al., “Experimental validation of a new ultrasoundmethod for the simultaneous assessment of radial and longitudinal
myocardial deformation independent of insonation angle,” Circulation, vol. 112, no. 14, pp. 2157–2162, Oct. 2005, doi: 10.1161/
CIRCULATIONAHA.105.554006.
[45] K. Saito et al., “Comprehensive evaluation of left ventricular strain
using speckle tracking echocardiography in normal adults: Comparisonof three-dimensional and two-dimensional approaches,” J. Amer. Soc.
Echocardiography, vol. 22, no. 9, pp. 1025–1030, Sep. 2009, doi:
10.1016/j.echo.2009.05.021.
[46] H. Hasegawa, “Phase-sensitive 2D motion estimators using frequency
spectra of ultrasonic echoes,” Appl. Sci., vol. 6, no. 7, pp. 1–19,
Apr. 2016, doi: 10.3390/app6070195.
[47] S. Nunome, R. Nagaoka, and H. Hasegawa, “Accuracy evaluation of 3D
velocity estimation by multi-frequency phase-sensitive motion estimator
under various specifications of matrix array probe,” Jpn. J. Appl. Phys.,vol. 59, no. SK, Jul. 2020, Art. no. SKKE01, doi: 10.35848/13474065/ab78e7.

Yu Obara was born in Iwate, Japan, in August
1996. He received the B.S. degree in electricalengineering and the M.S. degree in biomedical engineering from Tohoku University, Sendai,
Japan, in 2019 and 2021, respectively, where
he is currently pursuing the Ph.D. degree inbiomedical engineering.
His research interests include transcutaneous
measurements of the heart wall vibrations andmyocardial response to propagation of electricalpotential.
Mr. Obara is a member of the Acoustical Society of Japan and theJapan Society of Ultrasonics in Medicine.

1369

Shohei Mori (Member, IEEE) was born inMiyazaki, Japan, in February 1991. He received
the B.S., M.S., and Ph.D. degrees in engineeringfrom the Tokyo Institute of Technology, Tokyo,
Japan, in 2013, 2015, and 2017, respectively.
Since 2017, he has been an Assistant Professor with the Department of Electrical Engineering, Graduate School of Engineering, Tohoku
University, Sendai, Japan. His research interests include ultrasound measurements of organ
dynamics for cardiovascular diseases and quantification of the amplitude envelope statistics for tissue characterizationof liver diseases.
Dr. Mori is a member of the Acoustical Society of Japan and an
Engineering Fellow of the Japan Society of Ultrasonics in Medicine.

Mototaka Arakawa (Member, IEEE) was born in
Sendai, Japan, in January 1971. He received the
B.S., M.S., and Ph.D. degrees in electrical engineering from Tohoku University, Sendai, Japan,
in 1993, 1995, and 2000, respectively.
He is currently an Associate Professor with
the Graduate School of Biomedical Engineering, Tohoku University. His research interests
include viscoelasticity measurements of biological tissues and the development of acoustic
microscopy for tissue and cell characterization.
Dr. Arakawa is a member of the Acoustical Society of Japan, theJapan Society of Ultrasonics in Medicine, and the Institute of Electronics
Information and Communication Engineering of Japan.

Hiroshi Kanai (Senior Member, IEEE) wasborn in Matsumoto, Japan, in November 1958.
He received the B.E., M.E., and Ph.D. degrees
in electrical engineering from Tohoku University,
Sendai, Japan, in 1981, 1983, and 1986, respectively.
From 1986 to 1988, he was a Research Associate with the Education Center for Information
Processing, Tohoku University. From 1990 to
1992, he was a Lecturer with the Department
of Electrical Engineering, Faculty of Engineering,Tohoku University. From 1992 to 2001, he was an Associate Professor
with the Department of Electrical Engineering, Faculty of Engineering,Tohoku University. From 2012 to 2015, he was the Dean of the Graduate
School of Engineering, Tohoku University, where he was the Vice President of research consolidation and university reform from 2015 to 2018.
Since 2001, he has been a Professor with the Department of Electronic
Engineering, Graduate School of Engineering, Tohoku University, where
he has also been a Professor with the Department of Biomedical
Engineering, Graduate School of Biomedical Engineering, since 2008.
His present interests are in transcutaneous measurement of the heart
wall vibrations and myocardial response to propagation of electricalpotential and cross-sectional imaging of elasticity around atherosclerotic
plaque with transcutaneous ultrasound for tissue characterization of thearterial wall.
Dr. Kanai is a member of the Acoustical Society of Japan, a fellow of
the Institute of Electronics Information and Communication Engineeringof Japan, and a member of the Japan Society of Ultrasonics in Medicine,
the Japan Society of Medical Electronics and Biological Engineering, andthe Japanese Circulation Society. Since 1998, he has been a member of
the Technical Program Committee of the IEEE Ultrasonics Symposium.
Since 2011, he has been a Board Member of the International Congress
on Ultrasonics. Since 2013, he has been an Associate Editor of theIEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCYCONTROL. He was an Editor of the Journal of Medical Ultrasonics andthe Japanese Journal of Medical Ultrasonics from 2012 to 2018.

Authorized licensed use limited to: TOHOKU UNIVERSITY. Downloaded on April 06,2022 at 05:19:57 UTC from IEEE Xplore. Restrictions apply.