Japanese Journal of Applied Physics
Vol. 47, No. 5, 2008, pp. 4155–4164#2008 The Japan Society of Applied Physics

Automated Segmentation of Heart Wall Using Coherence Among Ultrasonic RF Echoes
Takaomi KINUGAWA, Hideyuki H ASEGAWA, and Hiroshi K ANAI
Graduate School of Engineering, Tohoku University, Sendai 980-8579, Japan(Received November 22, 2007; revised February 1, 2008; accepted February 12, 2008; published online May 23, 2008)

In conventional ultrasonic tomographic images, the heart wall cannot be distinguished from the cardiac lumen automaticallyon the basis of only the echogenicity. One of the biggest problems is that echogenicity, which corresponds to the amplitude ofan RF echo, in the heart wall is as low as that in the lumen. In this study, ultrasonic RF echoes from the heart wall and lumen
were analyzed in the frequency domain in order to distinguish the heart wall from the lumen automatically. Temporal changes
in complex frequency spectra were evaluated using the magnitude-squared coherence function. The coherence function of RF
signals scattered from the interventricular septum (IVS) was high. In contrast, the coherence function in the right ventricle(RV) and left ventricle (LV) was low because the scatterers (blood cells) slipped oﬀ from the focal area of the ultrasonic beamby blood ﬂow. For automated identiﬁcation of the heart wall using the coherence function, the optimal threshold T0 ð f Þ for thecoherence function should be determined. In this study, on the basis of the Bayes decision rule, the optimum value of T0 ð f Þ
was determined. The coherence function of the region near the anterior wall in the RV was as high as that in the IVS because
there are artifacts in the region near the anterior wall owing to echoes from the external tissue resulting from the sidelobe.
However, the artifacts can be reduced by removing the stationary component from RF echoes before evaluating the coherencefunction. In vivo experimental results show that the diﬀerentiation of the heart wall from the lumen was improved
signiﬁcantly using the proposed method. [DOI: 10.1143/JJAP.47.4155]
KEYWORDS: ultrasonic RF echo, heart, region identification, frequency analysis, coherence function

1.

Introduction

In recent years, the number of patients suﬀering fromcardiovascular diseases such as cardiac infarction has
increased. Cardiac catheterization, X-ray imaging, computerized tomography (CT), magnetic resonance imaging(MRI), and ultrasonography are employed for diagnosis ofthese diseases.1) Among these methods, ultrasonography issuperior to the others in terms of noninvasive and real-time
diagnosis. In addition to the morphological diagnosis of theheart, identiﬁcation of the heart wall is important for theevaluation of myocardial function2–4) because the region(= heart wall) to be analyzed must be determined.
In conventional ultrasonic tomographic images, the heart
wall cannot be distinguished from the cardiac lumen
automatically using only echogenicity. One of the biggest
problems is that echogenicity, which corresponds to theamplitude of the RF echo, in the heart wall is as low as thatin the lumen. Moreover, in the ultrasonic tomographic
images of the heart, there are many artifacts due to thereﬂection of sidelobes by the external tissue such as ribs.
Under these circumstances, in the conventional manual
segmentation of the heart walls in B-mode images, theregion of the heart wall is subjectively distinguished fromthe lumen on the basis of the a priori knowledge about theheart structure.
Nillesen et al. reduced the overlap between brightnesses
in the heart wall and lumen from 35.4 to 19.6% bysmoothing brightness in each local region of interest in
ultrasonic tomographic images.5) However, this method is
not eﬀective in the region with a low echo level inside theheart wall, and the heart wall is not completely distinguishedfrom the lumen.
On the other hand, the frequency characteristics of RF
echoes contain more information (amplitude and phase ateach frequency) on echoes from scatterers in comparison
with echo amplitudes as long as signal-to-noise ratio (S/N)


E-mail address: hkanai@ecei.tohoku.ac.jp

is not too small. Therefore, there is a possibility that thelumen and heart wall whose echogenicities are similar canbe diﬀerentiated using the frequency characteristics of RF
echoes from them. Lizzi et al. estimated the eﬀective
diameter of an ultrasonic scatterer on the basis of the slope
and intercept of the power spectrum normalized by that of
the echo from the ﬂat plate.6,7) Watson et al. and Waters
et al. classiﬁed arterial plaques in vitro on the basis of theslope of a normalized power spectrum.8,9) Hall et al.
evaluated the anisotropy of the frequency dependence ofbackscatters in the human myocardium.10) Georgia et al.
characterized breast tissue on the basis of time-frequency
analysis of RF signals.11,12) In this study, complex frequency
spectra of echoes from the heart wall and lumen wereanalyzed for the automated identiﬁcation of those regions.
2.

Experimental Methods

2.1 Temporal changes in complex spectra of RF echoes
In this study, we focus on the diﬀerence between temporalchanges in complex frequency spectra of RF echoes from theheart wall and the lumen. As shown in Fig. 1, in the heartwall, the phase changes of RF echoes during a few
milliseconds [1.69 ms (1 frame), 3.38 ms (2 frames), . . .]are very small when the motion of the myocardium istracked (compensated) because the myocardium stays in the
focal area of the ultrasonic beam, whereas small phasechanges occur owing to a change in the thickness of
myocardia. On the other hand, in the lumen, the motion of
scatterers cannot be compensated because blood cells slip oﬀ
from the focal area of the ultrasonic beam by blood ﬂow.
Therefore, phases of echoes would change even during a few
milliseconds. In this study, the magnitude-squared coherencefunction was employed for evaluating this diﬀerence invariance in phase changes of RF echoes.
2.2 RF data acquisition
As illustrated in Fig. 2(a), the interventricular septum(IVS), right ventricle (RV), and left ventricle (LV) were
scanned by an ultrasonic beam, which scans sparsely to4155

Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)

T. KINUGAWA et al.

phase change phase change phase change
lumen

myocardia
frequency [%]

heart
wall
blood cells
lumen

ultrasonic beamtime

a few [ms]: blood cells in the window− 180 0180180−180 0phase change during
a few milliseconds [deg]Fig. 1. Schematic diagram showing diﬀerence between phase changes of reﬂected waves during a few milliseconds in heart wall andlumen.

1
normalized amplitude

0.8

RV
IVS
LV
9 8

AV
LA

7
5 4 3 21
ultraso6nic
beams

Ao

0.6
0.4

1.10 μ s

0.2
0
−0.2
−0.4
−0.6
−0.8

(a)

− 20 dB

1.10 μ s 1540 m/s = 1.69 mm
0

1

2
3
time [μs]45Fig. 3. Waveform of input ultrasonic pulse measured using hydrophone.

the RF signal was 15 MHz and the frame rate in themeasurement was 592 Hz. Each RF data yn ðxÞ contains asmall bias signal owing to the ultrasonic diagnostic equipment. The bias should be removed because the direct current
is out of band. The bias is not constant along the directionof depth and is removed in every short period ðK þ 1Þx
as follows:
K=2
X
1y^n ðxÞ ¼ yn ðxÞ 
yn ðx þ ixÞ;ð2:1Þ
K þ 1 i¼K=2(b)Fig. 2. (a) Illustration of tomographic view of human heart. (b) B-mode
image of human heart in its longitudinal-axis view. Arrows in (a) show
ultrasonic beams and the region surrounded by the line in (a) and (b)shows the acquisition area of RF signals.

realize a high temporal resolution in a wide area.3,13) First,
RF data fyn ðxÞg along each of nine beams were acquired
using a 3.75 MHz sector-type probe of ultrasonic diagnostic
equipment (ALOKA SSD-6500), where n and x are frame(time) and depth, respectively. The sampling frequency ofwhere y^n ðxÞ is the RF signal after removing the regional bias
and x is the sampling interval (¼ 51:3 mm) in the directionof depth. The number K þ 1 of the sampled points used for
averaging is set to be 33 (1.69 mm) by referring to a pulse
length of insoniﬁed ultrasound (width at 20 dB of the
maximum of the envelope) of 1.69 mm in the soft tissue(sound velocity: 1540 m/s), as shown in Fig. 3.
2.3 Artifact reduction
RF data fyn ðxÞg contain the components of echoes fromthe external tissue such as ribs owing to sidelobes which

4156

Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)

T. KINUGAWA et al.

result in stationary artifacts. The artifact component was
reduced using the feature that the external tissue does not
move markedly with time as follows:
L=2
X
1
y~n ðxÞ ¼ y^n ðxÞ ð2:2Þ
y^nþi ðxÞ;L þ 1 i¼L=2where y~n ðxÞ is the RF signal after artifact reduction. The mean
absolute velocity of the myocardia in the IVS was 11 mm/s.
When myocardia in the IVS displace over one ultrasonicwavelength ( ¼ 513 mm) during L þ 1 frames, eq. (2.2)
does not aﬀect echoes from the heart wall markedly because
the phase of y^x ðnÞ rotates over 360 during L þ 1 frames andthe second term in the right hand side of eq. (2.2) for the
ﬂuctuating component is close to zero in comparison withthat for the stationary component. Therefore, L þ 1 should be
set at a value larger than 13.8 [corresponding to a period ofð=2Þ=11 mm/s = 23.3 ms] (condition 1). The wavelengthof  513 mm is derived from sound velocity (¼ 1540 m/s)
divided by the center frequency (3 MHz) of the insoniﬁed
pulse shown in Fig. 3 and the two-way propagation of theultrasound pulse between the probe and the object isconsidered. However, the artifact component also displaces
toward the surface of the body by a mean velocity of
1.2 mm/s in the data shown in Fig. 2. Therefore, L þ 1should be set as small as possible in order to reduce theartifact component eﬀectively (condition 2). Consequently,
L þ 1 was set to be 15 (25 ms) as the solution of the trade-oﬀ
problem between conditions 1 and 2. The mean displacement of the artifact component in L þ 1 frames is suppressed
to 30 mm (about 1/10 of the wavelength) when L þ 1 is 15.
Therefore, the artifact component can be reduced suﬃciently
and the reduction of the other components is suppressed.
2.4 Estimation of complex frequency spectrum
The complex frequency spectrum Yn;x ð f Þ in each frame n
and at each depth x was estimated by applying the short-time
Fourier transform to the RF signal y~n ðxÞ using a Hanningwindow. The window length was set at 1.64 mm (32 sampled
points) by referring to an insoniﬁed ultrasonic pulse lengthof 1.69 mm in the soft tissue (sound velocity: 1540 m/s), as
shown in Fig. 3. Then, the temporal changes of the spectra
were analyzed by applying the short-time Fourier transformto multiple frames. The window position was moved so as tofollow the motion of a point of interest using the phased
tracking method4,14) in order to realize the analysis of thesame scatterers during 100 frames, and the spectra fYn;x ð f Þg
in 100 frames were obtained for each position x. Using thephased tracking method, the motion of a point of interest is
traced without constraint of the sampling interval x along
the depth direction, and the displacement of the point
between two adjacent frames is obtained as a continuous
value. However, the window position should be assigned to adiscrete value depending on the sampling interval x. Thisdiﬀerence xn;x between the estimated position of the pointof interest and the discrete position of the window leads tothe undesirable phase shift n;x ð f Þ between the estimatedspectra, Yn;x ð f Þ and Ynþ1;x ð f Þ. This phase diﬀerence n;x ð f Þ
is given by
n;x ð f Þ ¼

4 fxn;x ;c0

ð2:3Þ

where c0 is the sound velocity of 1540 m/s. Consequently,
n;x ð f Þ can be compensated as follows:Y^n;x ð f Þ ¼ Yn;x ð f Þe jn;x ð f Þ ;ð2:4Þ

where Y^n;x ð f Þ is the phase-corrected complex frequency
spectrum.
2.5 Estimation of coherence function
Using the frequency spectrum Y^n;x ð f Þ in each frame, thepower of the time (frame) averaged frequency spectra
jEn ½Y^n;x ð f Þj2 was obtained. In addition, the magnitudesquared coherence function (MSCF) jx ð f Þj2 was evaluatedas follows:
jEn ½Y^2n;x
ð f ÞY^2nþN;x ð f Þj2
jx ð f Þj2 ¼; ðN ¼ 1; 2; . . .Þ;En ½jY^2n;x ð f Þj2 En ½jY^2nþN;x ð f Þj2 
ð2:5Þ
where En ½ and  are the time (frame) averaging during 100
frames and complex conjugate, respectively. The MSCF is afunction for evaluating the diﬀerence between the characteristics of phase changes of echoes from the heart wall and thelumen as described in §2.1. In the numerator of eq. (2.5), theabsolute value is taken after calculation of the mean ofcomplex cross-spectra. By this operation, the variance intemporal changes in the complex spectra between two RF
signals in (2n)-th and (2n þ N)-th frames during 100 frames
at each frequency is evaluated, i.e., the changes in phases ofcomplex spectra signiﬁcantly aﬀect this averaging operation.
The MSCF jx ð f Þj2 evaluates such change regardless of
the amplitude of RF echoes as long as S/N is not small
because it is normalized by the powers En ½jY^2n;x ð f Þj2  and
En ½jY^2nþN;x ð f Þj2 . In eq. (2.5), RF signals in (2n)-th and(2n þ N)-th frames were used to reduce computation time
because the MSCF between two RF signals in (2n)-th and(2n þ N)-th frames is nearly the same as that between two RF
signals in (2n þ 1)-th and (2n þ 1 þ N)-th frames because of
the employed high frame rate. Although the method of region
identiﬁcation using brightness, in which the phase of an echo
is not evaluated, was proposed,5) our proposed method using
complex spectra, which contain both amplitude and phase, ismore suitable for region identiﬁcation.
2.6 Determination of optimal threshold for MSCF
For the automated identiﬁcation of the heart wall using the
MSCF jx ð f Þj2 , the optimal threshold level T0 ð f Þ for the
MSCF was determined on the basis of the Bayes decision
rule.15,16) The a posteriori probability Pði j jx ð f Þj2 Þ, whichshows the probability that jx ð f Þj2 belongs to the class i , isshown as follows:
1Pð1 j jx ð f Þj2 Þ ? Pð2 j jx ð f Þj2 Þ ! jx ð f Þj2 2
;
2
ð2:6Þ
where the two classes 1 and 2 are the heart wall and thelumen, respectively. Using Bayes’ theorem, the a posteriori
probability Pði j jx ð f Þj2 Þ is shown as follows:Pði j jx ð f Þj2 Þ ¼

pðjx ð f Þj2 j i ÞPði Þpðjx ð f Þj2 Þ

ði ¼ 1; 2Þ; ð2:7Þ

where pðjx ð f Þj2 Þ, Pði Þ, and pðjx ð f Þj2 j i Þ are theprobability of occurrence of jx ð f Þj2 , a priori probability

4157

Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)

T. KINUGAWA et al.

t 0 −84 ms t 0 t 0+84 ms
(1)

(2)

5

(3)

10
depth [mm]

RV
IVS
LV
(a) M−mode image before artifact redcution

0
depth
10.26 mm

(α)
(β)
(γ)
(δ)
(ε)

15
20
25
30

(α)
RV

(β)
(γ)

IVS

(δ)
(ε)

LV

35
40
−0.05

tracking locus

0
amplitude

0510(b) M−mode image after artifact redcution(c) ECG

15

(α)
(β)
(γ)

RV

(δ)
(ε)

IVS

20
25
30

LV

35

0
(d) PCG00.05(e) RF signal
before artifact reduction

depth [mm]depth
10.26 mm

0

0

40
−0.05

1000
time [ms]5000amplitude

0.05(f) RF signal
after artifact reduction

Fig. 4. (Color online) (a) M-mode image of heart along beam 5 shown in Fig. 2(a) before artifact reduction. (b) M-mode image ofheart along beam 5 after artifact reduction. (c) Electrocardiogram (ECG). (d) Phonocardiogram (PCG). (e) RF signal y^n0 ðxÞ in rapid
ﬁlling phase before artifact reduction. (f) RF signal y~n0 ðxÞ in rapid ﬁlling phase after artifact reduction.

showing the occurrence of class i , and the conditional
probability density function showing the occurrence ofjx ð f Þj2 belonging to class i , respectively. By substituting
eq. (2.7) into eq. (2.6), the following eq. (2.8) is derived.
pðjx ð f Þj2 j 1 ÞPð1 Þ ? pðjx ð f Þj2 j 2 ÞPð2 Þ

1
:
ð2:8Þ
! jx ð f Þj2 2
2
In this study, the conditional probability density functionpðjx ð f Þj2 j i Þ and a priori probability Pði Þ were calculated on the basis of the histograms of jx ð f Þj2 of the heart
wall and lumen, which had been manually segmented. Usingeq. (2.8), we determined the optimal threshold T0 ð f Þ forjx ð f Þj2 .
3.

In vivo Experimental Results

3.1 RF data acquisition
RF data were acquired from the heart of a healthy 25year-old male, and the analysis described in the previous
section was applied.
Figure 2(b) shows a B-mode image. RF data wereobtained for the nine beams indicated by the arrows within
the acquisition area surrounded by the line in Fig. 2(a).
Figures 4(a) and 4(b) show M-mode images along beam 5
before and after artifact reduction, respectively. Figures 4(c)and 4(d) show the electrocardiogram (ECG) and phonocardiogram (PCG), respectively. The stationary artifacts in theregion () near the anterior wall shown in Fig. 4(a) werereduced, as shown in Fig. 4(b). The central time in theanalyzing period t0 was set so that it corresponds to one ofthree cardiac phases: (1) ejection, (2) rapid ﬁlling, and (3)reduced ﬁlling. RF data were analyzed around t0 from t0 
84 to t0 þ 84 ms. The RF signals y^n0 ðxÞ and y~n0 ðxÞ, which arerespectively before and after artifact reduction, at t0(¼ frame n0 ) in the rapid ﬁlling phase are shown in
Figs. 4(e) and 4(f), respectively. The amplitude of the RF
signal in the region () near the anterior wall decreased after
artifact reduction, as shown in Figs. 4(e) and 4(f). The
instantaneous position of each point of interest estimated by
the phased tracking method is indicated in Fig. 4(b) by theyellow lines of (), (), (), (), and (").
Regions (), (), (), (), and (") in Figs. 4(b) and 4(f) are
the positions of windows for RV (near the anterior wall),
RV, IVS (hyperechoic), IVS (hypoechoic), and LV, respectively. These ﬁve regions were chosen to compare thecharacteristics of echoes from the heart wall and lumen.
Note that the amplitudes of RF signals in regions (), (),and (") were similarly low. Figures 5(a) and 5(b) show thehistograms of brightness b^n ðxÞ and b~n ðxÞ before and after
artifact reduction, normalized by their maximum values, inthe heart wall and lumen acquired from nine beams shown in
Fig. 2(a), respectively. The brightness bn ðxÞ is deﬁned as thecommon logarithm of the envelope of the RF signal asfollows:
bn ðxÞ ¼ logfLPF½ðyn ðxÞÞ2 g;ð3:1Þ

where LPF½ shows the low-pass ﬁltering. The overlap ofbrightness b^n ðxÞ between the heart wall and lumen was

4158

probability density [%]probability density [%]Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)

30

T. KINUGAWA et al.

(a) before artifact reduction

15

lumen
overlap 40.2%
heart wall

0
0
30

0.2
0.4
0.6
0.8
normalized brightness

1(b) after artifact reduction

25
20

lumen

15
10
5
0
0

overlap 33.4%
heart wall
0.2
0.4
0.6
0.8
normalized brightness

1

Fig. 5. Histograms of brightness in heart wall and lumen acquired from 9
beams (a) before and (b) after artifact reduction. In this analysis, theregion near the anterior wall (<10:26 mm), in which the artifact was
dominant, was not included.

40.2% before artifact reduction. After artifact reduction, theoverlap of brightness b~n ðxÞ was 33.4%. These large overlaps
show that region identiﬁcation of the heart wall using thebrightness is diﬃcult.
3.2

Estimation of complex frequency spectrum andcoherence function
Figures 6(1-a) and 6(2-a) show the three instantaneous
power spectra jY^n;x ð f Þj2 and the power of the time-averaged
complex spectrum jEn ½Y^n;x ð f Þj2 along beam 5 before andafter artifact reduction, respectively, which were acquired
from the regions () to (") in the rapid ﬁlling phase shown in
Fig. 4(b). The shape of the power spectra changed with timeeven when the same region was analyzed. However, almost
all the power spectra and power of the temporally averaged
power spectra had the peaks at about 2 MHz ( fp ) and also
had the other peaks at 2 fp . In region () in RV, except for theregion near the anterior wall, and in region (") in LV, thepower of the averaged complex spectra was lower than the
instantaneous power spectra by 10 to 20 dB. Particularly inthe high-frequency range (¼ 2 fp ), the decrease in power of
the averaged complex spectra in the lumen was marked. In
the IVS [regions () and ()], however, the power of theaveraged spectra was similar to that of the instantaneous
spectra.
Figures 6(1-b) and 6(2-b) show the MSCF jx ð f Þj2
obtained during the period t0  84  t  t0 þ 84 ms in
the rapid ﬁlling phase at each depth x along beam 5 before
and after artifact reduction, respectively. In Fig. 6(b), theframe intervals N between two frames for the evaluation ofcorrelation are 1 (1.69 ms) and 2 (3.38 ms). The MSCF of the
RF echo scattered from the IVS was high. In contrast, theMSCF in the lumen was low because the scatterers (blood
cells) slipped oﬀ from the focal area of the ultrasonic beamcaused by blood ﬂow. As shown in Fig. 6(2-b), in the heartwall, the MSCF in the frequency range 2 fp was as high asthat in fp . On the other hand, in the lumen, the MSCF in thefrequency range 2 fp was lower than that in fp .
Figure 7(a) shows the original B-mode image. Figure 7(b)shows the MSCF jx ð f Þj2 in multiple frames in three cardiac
phases. In Fig. 7(b), the frame intervals N between twoframes for the evaluation of correlation are 1 (1.69 ms) and 2(3.38 ms). In every cardiac phase, the MSCF in the heart
wall was higher than that in the lumen.
The MSCF in the RV-lumen region () near the anteriorwall in Fig. 6(1-b) was as high as that in the IVS because theartifact, which did not move markedly with time, wasdominant in this region. In contrast, in Fig. 6(2-b), theMSCF in the region () near the anterior wall was reduced
using artifact reduction in eq. (2.2). However, the MSCF in
the heart wall did not change markedly between before andafter artifact reduction, as shown in Figs. 6(1-b) and 6(2-b).
Therefore, it was shown that eq. (2.2) was the appropriate
preprocessing for region identiﬁcation using the MSCF
because artifact reduction did not aﬀect the other echo
components signiﬁcantly.
3.3 Region identiﬁcation using MSCF
As shown in Fig. 5(b), the overlap of brightnesses b~n ðxÞ
between the heart wall and lumen was too large fordiﬀerentiation of the heart wall from the lumen using
brightness. On the other hand, as shown in Fig. 8, the overlap(%overlap) of the MSCF jx ð f Þj2 was only 3.0% in the rapid
ﬁlling phase at 4.2 MHz. These histograms correspond toprobability density functions fpðjx ðf Þj2 j i Þg ði ¼ 1; 2Þ andthe %overlap shows the probability of misclassiﬁcation whenusing eq. (2.8). In these analyses shown in Figs. 5 and 8, theregion () near the anterior wall (<10:26 mm), in which theartifact was dominant, was not included in order to determine
the optimal threshold T0 ð f Þ for the MSCF. Figure 9(a) showsthe %overlap of the MSCF between the heart wall and thelumen in three cardiac phases. The overlap of the MSCF at
the frequency 2 fp ( 4 MHz) was slightly smaller than that at
fp ( 2 MHz) and the overlap was the smallest when N ¼ 1.
Therefore, the frequency 2 fp and N ¼ 1 were employed in
the following region identiﬁcation.
On the basis of eq. (2.8), we determined the optimal
threshold T0 ð f Þ for the MSCF, as shown in Fig. 9(b).
Subsequently, the results of region identiﬁcation are shownin Fig. 10, where the regions shown in red are identiﬁed as
the heart wall and those shown in blue are identiﬁed as thelumen. In this ﬁgure, white lines show the manually
determined boundary between IVS and ventricles. Theselines are determined by referring to the leading edges of RF
echoes from interfaces. Although there were misclassiﬁcations around the boundary between IVS and ventricles
because the analyzed window has a width of 1.64 mm, the
IVS is clearly diﬀerentiated from the RV and the LV in all
cardiac phases, and the opening and closing of the aortic
valve (AV) can also be observed in Fig. 10.
The same procedure was applied to another subject

4159

Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)

25

(α)

(b)

50
25

0

0

−25

−25

power spectrum [dB]−50−50500 1 2 3 4 5 6 7 500 1 2 3 4 5 6 7
(β)
(ε)
25
25
0

0

−25

−25

−50
−50
0 1 2 3 4 5 6 7500 1 2 3 4 5 6 7
frequency [MHz]
(γ)
25
0
−25
−50
0 1 2 3 4 5 6 7
frequency [MHz]power spectrum [dB]

(a)

50

(α)

N =1

(δ)

N =2
(α)

central depth of windowed RF signal [mm]

50

10

10

50

25

25

0

0

−25

−25

(γ)

20

0

0

−25

−25

−50 0 1 2 3 4 5 6 7 −500 1 2 3 4 5 6 750
frequency [MHz]
(γ)
25
0

20
IVS

(δ)
30

30
(ε)

LV

0 1 2 3 4 5 6 7
frequency [MHz](1) before artifact reduction
(b)
N =1
(δ)
(α)

−50 0 1 2 3 4 5 6 7 −500 1 2 3 4 5 6 750(ε)(β) 50
25
25

RV

(β)

40

central depth of windowed RF signal [mm](a)T. KINUGAWA et al.

1040
0 1 2 3 4 5 6 7
frequency [MHz]
N =2

RV

10
(β)
(γ)

20

20
IVS
(δ)

30

30
(ε)
LV

−25
−50 0 1 2 3 4 5 6 7
frequency [MHz]4040
0 1 2 3 4 5 6 7
0 1 2 3 4 5 6 7
frequency [MHz]frequency [MHz](2) after artifact reduction
magnitude−squared
instantaneous power spectrum
coherence function
t0 +50 [ms]t0 −50 [ms]t0 [ms]0

power of time−averaged comprex spectrum

0.5

1Fig. 6. (Color online) (a) Power spectra jY^n;x ð f Þj2 of RF signal along beam 5 and power of time-averaged complex spectrum
jEn ½Y^n;x ð f Þj2 . (b) MSCF jx ð f Þj2 between spectra in 2n-th and (2n þ N)-th frames along beam 5. (1) Before artifact reduction. (2)After artifact reduction.

(healthy 24-year-old male) in the rapid ﬁlling phase.
Figure 11(a) shows the original B-mode image. The brightness inside IVS was very low, the artifact in the RV was large
in this heart, and the overlap of the brightness b~n ðxÞ between
the heart wall and lumen was 47.0%. Figure 11(b) shows animage of the MSCF (4.2 MHz, N ¼ 1). Figure 11(c) shows
the result of region identiﬁcation. Although there was
misclassiﬁcation around the aorta (Ao), the IVS was clearly
diﬀerentiated from the RV and the LV.

4.

Discussion

4.1 Relationship between phase variance and MSCF
The possibility of misclassiﬁcation is the smallest at
frequency 2 fp as described in §3.3. In this section, the reason
is discussed using the relationship between the variance inphase variations and the MSCF.
In the analyzed period (100 frames), a larger variance2
ð f Þ of the phase variation  n ð f Þ between the spectra

4160

Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)

T. KINUGAWA et al.

(1) ejection phase10.26 mm
IVS

N=1

RV
LV

Ao

(2) rapid filling phase

N=2

N=1

N=2

(3) reduced filling phase

N=1

N=2
2.3
(a) original B−mode

3.5
frequency [MHz](b) Images of MSCF

0

4.7
0.5

1

probability density [%]Fig. 7. (Color online) (a) Original B-mode images. (b) Images of MSCF jx ð f Þj2 at each frequency in 3 cardiac phases. Color codes are
the same as those in Fig. 6(b).

8heart wall
4 lumenoverlap 3.0%
0
0
0.2
0.4
0.6
0.8
1
magnitude− squared coherence function

Fig. 8. Histograms of MSCF jx ð f Þj2 (frequency: 4.2 MHz, N ¼ 1) inheart wall and lumen acquired from 9 beams. In this analysis, the region
near the anterior wall (<10:26 mm), in which the artifact was dominant,
was not included.

Y^2n;x ð f Þ and Y^2nþN;x ð f Þ in (2n)-th and (2n þ N)-th frames
leads to a lower MSCF. The phase variation  n ð f Þ isshown as follows: n ð f Þ ¼ arg½Y^2n;x
ð f ÞY^2nþN;x ð f Þ;

ð4:1Þ

The relationship between the standard deviation of thepﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2 ð f Þ of the phase variation  ð f Þ between the spectranY^2n;x ð f Þ and Y^2nþN;x ð f Þ and MSCF jx ð f Þj2 was simulated

under the assumption that the power spectrum does notvary with time but only the phase varies. The distributionof  n ð f Þ was assumed to be Gaussian whose average
is zero deg, as shown in Fig. 12(a). The number of½Y^2n;xð f ÞY^2nþN;x ð f Þ was set at 10,000, which is suﬃcientfor obtaining smooth graphs. Figure 12(b) shows thepﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
result2ð f Þ
of simulation to estimate the relationship between
2and jx ð f Þj . When it is assumed that the phase of the RF
signal depends on the location
 n ð2 fp Þ ispﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃof scatterers,
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2 ð2 f Þ is twice
2 ð f Þ because
twice  n ð fp Þ and alsopp
 n ð f Þ occurs owing
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃto the displacement of scatterers.
2 ð f Þ is as small as 10 , as is expected
Therefore, whenpin the heart wall, jx ð fpp
Þj2ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
is 0.95 and jx ð2 fp Þj2 is 0.8. On
2the other hand, whenð fp Þ is 30 , as is expected in2
the lumen, jx ð fp Þj is 0.6 and jx ð2 fp Þj2 is signiﬁcantly
decreased to 0.15. Therefore, the ratio between the MSCF in
the heart wall and the lumen at 2 fp is larger than that at fp .
Figure 13 shows the distribution of the phase variation
 n ð f Þ in in vivo data obtained from regions (), (), (),and (") in the rapid ﬁlling phase shown in Fig. 4(b).
Figures 13(a) and 13(b) show the results at 2.3 MHz ( fp )and 4.5 MHz ( 2 fp ), respectively. From these results, theassumption that the distribution of  n ð f Þ is Gaussian
whose average is zero deg is shown to be adequate.
Therefore, the complex coherence function x ð f Þ was not
employed but the MSCF jx ð f Þj2 was employed in this study

4161

Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)

T. KINUGAWA et al.

(a)
RV
IVS

overlap [%]100 (a)cardiac phase90(1) (2) (3)
80
N =1
N =2
70
N =3
60
50
40
30
20
10
00 1 2 3 4 5 6
frequency [MHz]

Ao
LV

AV

(b)
1

7

0.5
0
(c)

0.8

RV

0.6

IVS

0.4
LV

0.2
1

2 3 4 5 6
frequency [MHz]7

Fig. 9. (a) Overlaps of MSCF jx ð f Þj2 between heart wall and lumen. (b)Optimal threshold levels T0 ð f Þ of MSCF jx ð f Þj2 determined by MSCFs
in manually assigned heart wall and lumen. (1) Ejection phase. (2) Rapid
ﬁlling phase. (3) Reduced ﬁlling phase.

Fig. 11. (Color online) Results for 24-year-old male. (a) Original B-mode
image. (b) Image of MSCF in rapid ﬁlling phase (4.2 MHz, N ¼ 1). (c)Region identiﬁcation using optimal threshold [T0 ð f Þ ¼ 0:47, f ¼
4:2 MHz].

3(a)

T0 ( f )=0.79( f =4.2 MHz)

1.5

0
0
60 120 180−180 −120 −60
phase variation Δφn ( f ) [deg]
1
(b)
0.9
fp
0.8
0.7
2 fp
0.6
0.5
fp
0.4
0.3
0.2
0.1
2 fp
0
0 10 20 30 40 50 60 70 80 90 100
σφ2( f ) [deg]

(b)

AV

T0 ( f )=0.68( f =4.2 MHz)(c)

T0 ( f )=0.61( f =4.2 MHz)
Fig. 10. (Color online) Region identiﬁcation using optimal thresholddetermined in Fig. 9(b). (a) Ejection phase. (b) Rapid ﬁlling phase. (c)Reduced ﬁlling phase.

because the phase of x ð f Þ is 0 when the average of  n ð f Þ
is zero deg and echo amplitude is constant along thedirection of frame.
4.2 Relationship between N and MSCF
The possibility of misclassiﬁcation is the smallest when
N ¼ 1 as described in §3.3. In this section, the relationshipbetween N and the MSCF is discussed using the relationship(a)probability density [%]

00

Ao
AV

magnitude−squared coherence function

threshold level1 (b)σ ( f )=10°2
φσφ2( f )=20°
σφ2( f )=40°
σφ2( f )=80°

Fig. 12. (a) Gaussian distribution. (b) Simulation of relationship between
magnitude-squared
coherence function jx ð f Þj2 and standard deviation
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
p
2 ð f Þ of phase variation  ð f Þ between spectrum Y^n2n;x ð f Þ in 2n-th
frame and Y^2nþN;x ð f Þ in (2n þ N)-th frame.

between blood ﬂow velocity and the MSCF and that between
the change in thickness of the myocardia and the MSCF.
The maximum blood ﬂow velocity in the LV measured
by the ultrasonic Doppler method was about 0.7 m/s inthe rapid ﬁlling phase [Fig. 14(a)]. The ultrasonic beam

4162

Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)(a) 2.3 MHz80(β) RV lumen

60

60

40

40

20

20

0
−180 −120 −60
80

0

60

0
120 180 −180 −120 −6080frequency [%](γ) IVS (hyperechoic)

60

40

40

20

20
0

60

0
120 180 −180 −120 −60
80

60

60

40

40

20

20
0

60

120 1800

60

120 1800120 180 −180 −120 −60
80

0

60

120 180

60

60

40

40

20

20
0

60

0
120 180 −180 −120 −60

phase variation Δφn( f ) [deg]060

120 180

Fig. 13. Distribution of phase variation  n ð f Þ between spectra in (2n)-thand (2n þ N)-th frames (N ¼ 1). (a) 2.3 MHz. (b) 4.5 MHz.

blood flow velocity(a)amplitude(b)+

0.7 m/s (rapid filling phase)+0.7 m/s
+0.5
0

−

100
90
80
70
60
50
40
30
20
10
0

0

20
40
60
80 100 120 140
frame number of time averaging

Fig. 15. Relationship between %overlap of MSCF and frame number oftime averaging.

(ε) LV lumen(ε) LV lumen

0−180 −120 −6060(δ) IVS (hypoechoic)(δ) IVS (hypoechoic)

0−180 −120 −60800(γ) IVS (hyperechoic)600−180 −120 −6080(b) 4.5 MHz(β) RV lumen
overlap [%]80

T. KINUGAWA et al.

−0.5

PCG
ECG

1
0.9
0.8
0.7
0.6
1.45 mm
0.5
0.4
0.3
0.2
0.1
0
− 2 − 1.5 − 1 − 0.5 0 0.5 1 1.5 2
displacement along radial directionof beam [mm]Fig. 14. (a) Blood ﬂow velocity in LV measured by ultrasonic Doppler
method. (b) Maximum of amplitude of echo from ﬁne wire. Ultrasonicbeam width was determined by width at half maximum.

diameter measured as a width at half maximum of theamplitude of an echo from a ﬁne metal wire was 1.45 mm[Fig. 14(b)]. Therefore, blood cells may completely slip oﬀ
from the beam during 2 ms (1.2 frames) in the rapid ﬁllingphase. That is why, in the lumen, the MSCF is slightly highwhen N ¼ 1 but is almost zero when N ¼ 2 in the cardiac
phases in which blood ﬂow is fast as in the ejection and rapid
ﬁlling phases. In the reduced ﬁlling phase, the MSCF in thelumen is also slightly high when N ¼ 2 because blood ﬂow
is slow (<0:3 m/s) and the scatterers stay in the same focal
area of the beam during a few frames. On the other hand, in
the IVS, the MSCF also decreases in the ejection and rapid
ﬁlling phases when N is increased because the number andthickness of myocardial ﬁbers included in an analyzed
window varies during N frames. In the reduced ﬁllingphase, the MSCF in the heart wall does not decrease
markedly when N is increased because the number andthickness of myocardium ﬁbers included in a analyzing
window does not vary markedly during N frames. Therefore,
the optimal threshold levels of the MSCF are almost thesame in three cardiac phases when N ¼ 1 but those in
the ejection and rapid ﬁlling phases are much lower than thatin the reduced ﬁlling phase when N ¼ 2, as shown in
Fig. 9(b). Moreover, in Fig. 9(a), %overlaps of the MSCFs
in the heart wall and lumen at N ¼ 1 are smaller than thoseat N ¼ 2 in the ejection and rapid ﬁlling phases because thedecrease in the MSCF owing to the increase in N in the heart
wall is larger than that in the lumen. Therefore, in this study,
the optimal value of N was determined to be N ¼ 1.
4.3 Number of frames in time averaging
In eq. (2.5), the number of frames in time averaging was
set to be 100 frames. This setting is shown to be appropriate
as described below.
Figure 15 shows the relationship between the %overlap of
the MSCF and the number of frames in time averaging. A
small number of frames in time averaging would take a large
MSCF also in the lumen. On the other hand, a large numberof frames in time averaging would take a small MSCF in theheart wall because of the large change in the echo signal
2
ð f Þ owing to a large change in the thickness of the
myocardia during the averaging period. When the number offrames in time averaging is about 100, %overlap becomes
the smallest, as shown in Fig. 15. Therefore, in this study,
100-frame averaging was employed.

4163

Jpn. J. Appl. Phys., Vol. 47, No. 5 (2008)

5.

T. KINUGAWA et al.

Conclusions

On the basis of the evaluation of the temporal changes in
the frequency spectra of RF echoes from the heart using themagnitude-squared coherence function, the diﬀerence between time changes of echoes from the heart wall and lumen
was shown. Diﬀerentiation of the heart wall from the lumen
was improved by artifact reduction and use of the highfrequency components (at 2 fp ). The possibility of an
automated region identiﬁcation of the heart wall using themagnitude-squared coherence function was shown by thein vivo experimental results that the heart wall can be clearly
diﬀerentiated from the lumen when the optimal threshold forthe magnitude-squared coherence function is determined.

1) F. H. Netter: The Ciba Collection of Medical Illustrations Heart (CibaGeigy Japan, Takarazuka, 1975) pp. 22 – 47.
2) H. Kanai, H. Hasegawa, and K. Imamura: Jpn. J. Appl. Phys. 45(2006) 4718.
3) H. Yoshiara, H. Hasegawa, H. Kanai, and M. Tanaka: Jpn. J. Appl.
Phys. 46 (2007) 4889.
4) H. Kanai, H. Hasegawa, N. Chubachi, Y. Koiwa, and M. Tanaka:
IEEE Trans. Ultrason. Ferroelectr. Freq. Control 44 (1997) 752.

5) M. M. Nillesen, R. G. P. Lopata, I. H. Gerrits, L. Kapusta, H. J.
Huisman, J. M. Thijssen, and C. L. de Korte: Ultrasound Med. Biol. 33(2007) 1453.
6) F. L. Lizzi, M. Greenebaum, E. J. Feleppa, and M. Elbaum: J. Acoust.
Soc. Am. 73 (1983) 1366.
7) F. L. Lizzi, M. Ostromogilsky, E. J. Feleppa, M. C. Rorke, and M. M.
Yaremko: IEEE Trans. Ultrason. Ferroelectr. Freq. Control 34 (1987)
319.
8) R. J. Watson, C. C. McLean, M. P. Moore, T. Spencer, D. M. Salter, T.
Anderson, K. A. A. Fox, and W. N. McDicken: Ultrasound Med. Biol.
26 (2000) 73.
9) K. R. Waters, S. L. Bridal, C. Cohen-Bacrie, C. Levrier, P. Fornès, and
P. Laugier: Ultrasound Med. Biol. 29 (2003) 1521.
10) C. S. Hall, E. D. Verdonk, S. A. Wickline, J. E. Perez, and J. G. Miller:
J. Acoust. Soc. Am. 101 (1997) 563.
11) G. Georgia and F. S. Cohen: IEEE Trans. Ultrason. Ferroelectr. Freq.
Control 48 (2001) 355.
12) G. Georgia, F. S. Cohen, C. W. Piccoli, F. Forsberg, and B. B.
Goldberg: IEEE Trans. Ultrason. Ferroelectr. Freq. Control 48 (2001)
364.
13) H. Kanai: IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52 (2005)
1931.
14) H. Kanai, M. Sato, Y. Koiwa, and N. Chubachi: IEEE Trans. Ultrason.
Ferroelectr. Freq. Control 43 (1996) 791.
15) H. Kanai: Spectrum Analysis of Sound and Vibration (Corona, Tokyo,
1999) p. 97.
16) R. J. Schalkoﬀ: Pattern Recognition: Stastical, Structual and Neural
Approaches (John Wiley & Sons, New York, 1992) p. 298.

4164