IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

791

Transcutaneous Measurement and Spectrum
Analysis of Heart Wall Vibrations
Hiroshi Kanai, Member, IEEE, Michie Sato, Yoshiro Koiwa, and Noriyoshi Chubachi, Member, IEEE

Abstract—For the noninvasive diagnosis of heart disease basedon the acoustic and elastic characteristics of the heart muscle, itis necessary to transcutaneously measure small vibration signals,
including components with an amplitude of less than 100 m,from various parts of the heart wall continuously for periods ofmore than several heartbeats in a wide frequency range up to 1
kHz. Such measurement, however, has not been realized by any
ultrasonic diagnostic methods or systems to date. By introducing
the constraint least-square approach, this paper proposes a new
method for accurately tracking the movement of the heart wall
based on both the phase and magnitude of the demodulated
signal to determine the instantaneous position of the object sothat the vibration velocity of the moving object can be accurately
estimated. By this method, small vibrations of the heart wall withsmall amplitudes less than 100 m on the motion resulting froma heartbeat with large amplitude of 10 mm can be successfully
detected with sufficient reproducibility in the frequency range
up to several hundred Hertz continuously for periods of about
10 heartbeats. The resultant small vibration is analyzed not
only in the time domain, but also in the frequency domain. As
confirmed by the preliminary experiments herein reported, the
new method offers potential for research in acoustical diagnosisof heart disease.

I. INTRODUCTIONTO diagnose the cardiovascular system based on the acoustic characteristics of the myocardium or the arterial walls,
it is necessary to measure small vibration signals on thesewalls in a wide frequency range and to accurately analyze
the spectrum of the resultant waveform. For the latter purposeof spectrum analysis, we have already developed a method
for estimating the spectral transition of heart sounds between
short-length signals of succeeding frames in low SNR cases[1]. By applying this method to the analysis of multiframe
signals of the fourth heart sound obtained during a stress test,
significant differences between the spectral transition patternsof 10 patients with myocardial infarction and 10 normal
subjects were clearly detected. However, when measuring aheart sound from the chest wall of a patient, it is difficult to
avoid the influence of the transfer characteristics of the lung.
Thus, the heart sounds were detected by an accelerometer
placed in the esophagus [1]. However, esophageal probe

Manuscript received March 14, 1995; revised February 28, 1996. This work
was supported in part by a Grant-in-Aid for Scientific Research (06750430,
06555113) from the Ministry of Education, Science and Culture of Japan andin part by Torei Science Foundation and Ono Science Foundation.
H. Kanai, M. Sato, and N. Chubachi are with the Department of ElectricalEngineering, Faculty of Engineering, Tohoku University, Sendai 980-77 Japan(e-mail: hkanai@chubachi.ecei.tohoku.ac.jp).
Y. Koiwa is with the First Department of Internal Medicine, School of
Medicine, Tohoku University, Sendai 980-77 Japan.
Publisher Item Identifier S 0885-3010(96)06340-X.

placement is not well tolerated in unanesthetized patients.
Moreover, it is difficult to detect local vibrations of small areasof various parts of the ventricle wall with such placement. It
is, therefore, desirable to use a probe placed on the chest wall
to transcutaneously detect small vibrations for each local areaby means of ultrasound.
Since motion detection based on the continuous-wave ultrasonic Doppler method was first demonstrated [2], numerous
elaborate techniques have been proposed for acoustic measurement of the velocity of the blood flow in the heart or thearteries based on the Doppler effect [3]–[11].
However, to measure heart wall vibration signals with small
amplitudes (less than 100 m) up to at least a few hundred
Hertz, the vibrations are superimposed on the motion withlarge amplitudes (about 10 mm) of several Hertz resulting
from the heartbeat during one cardiac cycle. Thus, there arelarge fluctuations (about 8 s) in the transit period of an
ultrasonic wave traveling from a transducer to the heart walland back during the period of one heartbeat. It is, therefore,
necessary to accurately identify the instantaneous position of
the heart wall.
In this paper, the term “vibration” is employed to describe
the rapid motion of the heart wall with small amplitude
including high-frequency components, which is difficult torecognize by M (motion)-mode echocardiography. Thus, avibration is distinguished from a motion with large amplitude
and low frequency, the displacement of which can be displayedin the M-mode image. A small vibration on the heart wall,
which is measured up to several hundred Hertz in this paper, isanalyzed not only in the time domain, but also in the frequency
domain, whereas the motion of the tissue or organs depicted instandard echocardiography is analyzed only in the time domain
because it has only low-frequency components.
Various ultrasonic techniques for measurement of the motionof tissue or organs using ultrasound have been reported in theliterature. The most widely employed technique is M-mode
echocardiography, in which the amplitude of each ultrasonicpulse shows the slow motion of the structure along theultrasonic beam. However, the spatial resolution is limited toseveral wavelengths, namely, more than 1 mm for ultrasoundof 3 MHz. The limit is determined by the length of theultrasonic RF pulse radiated by the ultrasonic transducer.
Moreover, to determine the displacement of organs as waveforms to be analyzed, it is necessary to trace the resultant wallmotion. However, this is not easily realized.
Several methods have been developed to measure changes
in the diameter of the arterial walls by tracking arterial

0885–3010/96$05.00  1996 IEEE

792

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

wall motion, including only low-frequency components withsmall amplitudes for detection of blood flow [12], [13], forevaluation of elastic properties of the arterial walls [14]–[16],for measurement of pressures [17], and for monitoring fetal
breathing [18].
A method has also been proposed and applied to measurearterial wall motion having large amplitudes in which the zerocrossing point of the echo signal from an arterial wall istracked for each transmitted RF pulse using the phase-locked
loop (PLL) technique [19]–[26]. In the method using the PLL
circuit, however, only the phase of the received RF signal
is considered and it is easily affected by additive noise as
pointed out in [27]. With each of these devices using the
PLL technique, experiments have shown that it is difficult
to lock onto and remain locked to the desired echo [27].
However, the reason has not been considered theoretically.
As will be pointed out in Section II-D of this paper, thismethod has an essential drawback in that the displacement
is coarsely digitized at the sampling rate determined by theclock frequency employed in the PLL circuit. An alternative
open-loop approach has been developed recently in which
the phase is detected at a fixed distance from the transducer to trace the motion of the layers [27]–[29]. In thesemethods, an ultrasonic transducer is directly attached to theepicardium.
Since the amplitude of the motion of the vessel walls is not
so large, the walls do not deviate from the selected sample
volumes during the cardiac cycle. For an artery with low
elasticity, however, the vessel walls deviate from the selected
sample volume. This disadvantage has been overcome byletting the sample volume track the vessel walls based on theassessed displacement [30]. The displacement waveforms andthe distension waveforms of the arterial walls are obtained andthen analyzed in the time domain. However, transcutaneous
measurement of the velocity signal on the heart wall has not
yet been realized. Moreover, the obtained waveforms have
not yet been analyzed as to their spectrum in the frequency
domain.
We have proposed, therefore, a noninvasive transcutaneous
method for measuring small vibrations on a local area of theheart wall [32]. The demodulated signal of the received signal
is A/D converted at a sampling frequency of 1 MHz. Fromthe phase of the received echoes, a small vibration signal ismeasured noninvasively using the ultrasonic Doppler method.
In this method, however, the object position is determined only
from the local maximum of the magnitude in the demodulated
received echo, and the spatial resolution in the tracking of theobject position highly depends on the sampling period of thedemodulated received signal.
This paper proposes a new method for accurately tracking
the heart wall movement using both the phase and magnitudeof the digitized demodulated signal to determine the instantaneous position of the heart wall, after which the velocity signal
is estimated. The resultant small vibration is analyzed not onlyin the time domain, but also in the frequency domain. Afterevaluating the validity of the proposed method by experiments
using a water tank, the proposed method is applied to the invivo detection of small vibration signals on the heart walls in

Fig. 1. Schematic representation of the transcutaneous measurement of aheart wall vibration, v t , from the chest wall using pulsive ultrasound basedon the standard procedure when the object position does not vary so much.
f0 are transmitted at a time
RF pulses with angular-frequency !0interval of T from an ultrasonic transducer on the chest surface. The
ultrasonic pulse reflected by the heart wall is received by the same ultrasonictransducer. The output signal, z t , is amplified and quadrature demodulation
is applied to the signal. The resultant complex signal is A/D converted ata sampling period of Ts and then separated into response signals fy x t g
for each transmitted pulse at a time t, where x denotes the distance of theobject from the ultrasonic transducer and where the spatial resolution xs inthe depth direction is equal to Ts 1 c0 = . Since the phase  x t of theresultant sectional complex signal y x t is given by !0 x t =c0 , fromthe phase difference  x t between the demodulated signals y x t andy x t
T of the received signals for the successively transmitted pulses
in the interval T , the average velocity v t of the object during the period
T is given.

()1( ; +1 )
1
1

=2

()

1(; )

2
(; )

()

(; )
1
(; )
2 ()
(; )

humans in order to show the possibility of acoustic diagnosisof the myocardium.
II. PRINCIPLE OF MEASUREMENT OF
SMALL VIBRATIONS ON A MOVING OBJECT
A. Standard Velocity Measurement
are transmitRF pulses with angular-frequency
ted at a time interval of
from an ultrasonic transducer
on the chest surface as shown in Fig. 1. Defining the acoustic
velocity as , the instantaneous distance of a moving object,
from the ultrasonic transducer is denoted bywhere
is the instantaneous period required for one-way
transmission from the ultrasonic transducer to the object. The
ultrasonic pulse reflected by the object is received by the same
, is amplified andultrasonic transducer. The output signal,
quadrature demodulation is applied to the signal. The resultant
complex signal is A/D converted at a sampling period ofand then separated into the response signals
for eachtransmitted pulse at a time , where denotes the distanceof the object from the ultrasonic transducer, and the spatialin the depth direction is equal toasresolution
shown in Fig. 1.

KANAI et al.: HEART WALL VIBRATIONS

793

The phaseof the resultant sectional complex signal
is given by the angular frequency
multiplied byin the one-way propagation from thetwice the delay time
ultrasonic transducer to the object as follows:

At first, based on the standard least-square method, let
us define the normalized mean squared difference between
and thethe demodulated complex signal
bypreceding signal(1)
The phase difference
between the demodulated signals
andof the received signals for theis given by
successively transmitted pulses in the interval(2)is the movement of thewhere
at a time , where it can be assumed
object in the periodthat the received interval almost coincides with the transmittedof the ultrasonic RF pulses. By dividing theinterval
movementby the period, the average velocity, denoted, of the object during the period
is givenby
between the successively
by the phase differenceandas follows:
demodulated signals(4)is the movement, to be determined, of the object
whereof the succeeding pulses, is the range
during the periodaround the previous object-position
where
is the change
the above difference is evaluated, andto.
in complex amplitude fromcorresponds to the length of the RF pulse
The rangetransmitted from the ultrasonic transducer. From (A9) inis achieved by
Appendix A, the minimization of

(5)

where

denotes the complex conjugate. For this value of, the minimum
is given by(3)
This standard procedure can be employed only when the object
does not vary much.
position
B. Tracking Object Position Using the Standard
Cross-Correlation Function
Since the object positionchanges by more than 10 mm,
for example, for the heart wall during one heartbeat period,
it is necessary to accurately track the instantaneous positionof the object for estimation of the vibration signalof the object. The transmitted ultrasonic pulse has a finite
length of several microseconds when the ultrasonic frequency
is several megahertz. Thus, there is a pulse train consisting ofa few discrete points at which each reflective signal has large
magnitude as shown in the right-hand side of Fig. 1. Therefore,
from theit is not easy to determine the object position. Moreover, there
magnitude of the demodulated signal
is a large phase difference between the adjacent points ina train when the object is moving. Thus, an error in thedetermination of the object position introduces a large errorin (3). Therefore, it
into the resultant vibration signal
is essential to accurately determine the instantaneous object.
positionIn this paper, in order to accurately determine both phaseduring the period
and the instantaneousshift, the following complex correlation is
object position
introduced in the determination procedure of the movement,, of the object position
from the preceding sectional.
data during the period(6)
from (A11) in Appendix A. The termin the numerator of (6) corresponds to themagnitude of the complex cross-correlation function between
and
.
and
However, when the received complex signals
are described by the damped sinusoidal signals asshown in Fig. 2, the above standard correlation procedure doesof the signal
not uniquely determine the delay
from
as described below. Let us assume that theandare digitizedreceived complex signals
and are given by the following simple model:

(7)
(8)

794

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

Fig. 2. The waveforms of a simple
and y x t
complex signals y x t
imaginary components are indicated by

(; )

model

of

the

received( ; + 1T ). The real andand 2, respectively. (a)= y(x; t)j = = z0 1 w(n) = e00 1 fcos (2n=9) +y1 (n) def
j sin (2n=9)g 1 w(n) in (7). (b) y2 (n) def= y(x; t + 1T )j = =z0 + 1w(n+n0 ) = e00 1( +5) fcos[2(n+5)=9]+j sin [2(n+5)=9]g1
w(n + 5) in (8).
x

n

: n

x

n

where
from

n

n

:

n

n

is the complex constant,is the delay of, andis the step function defined by
(9)

is
and
, the waveforms ofandare shown in Fig. 2(a) and (b), respectively.
For this simple model, the minimumin (6) of thenormalized mean squared difference between these signals is, the real constant ,calculated using the real function
and the complex function
derived in Appendix B. Fromof (6) is given by(A10), the minimum
When

( )

Fig. 3. (a) The characteristics of the minimum n MIN x of (11) and (12),which are derived from the standard normalized mean squared differencen0,
between the signals in Fig. 2(a) and (b). For the case x
is zero. Thus, it is impossible to uniquely determine the truen MIN x
value n0 of x from n MIN x . (b) The characteristics of the minimum
MIN x of (19) and (20), which are derived by introducing the restriction
into the normalized mean squared difference between the signals in Fig. 2(a)and (b). The MIN x takes the minimum only at the true lag of. Thus,the introduction of the restriction is effective in determination of lag between
these complex signals.

( )
0
( )

( )

( )

 0 = 05
05

C. Tracking Object Position Using theConstraint Least-Squares Method
As described above for the simple example, the waveformfits well withfor many candidates ofif the fitness is evaluated by the normalized mean squared
difference defined in (4). To solve the above problem, let us
introduce the constraint least-squares method [31] by which
the phase change rather than the magnitude change between
these signals is evaluated. Let us redefine the normalized meansquared difference of (4) between the demodulated complex
and the preceding signal
bysignalof(10), by substituting (B2), (B4), and (B6)
For the casein Appendix B into (10)

(13)
(11)
For the caseinto (10), by substituting (B3), (B4), and (B7)(12)
of (11) and
The characteristics of the minimum, by(12) are shown in Fig. 3(a). For the casein (4) by,is always zero.
settingThus, it is difficult to uniquely determine the true valueoffrom the curve ofobtained by the standard
complex correlation approach.

whereis the phase change fromto
. The magnitude ofis restricted to be one. This
restriction is realized by introducing the Lagrange multiplier
into the second term in the right-hand side of (13). As
under
described in Appendix C, the minimization ofis achieved by
the restriction(14)

KANAI et al.: HEART WALL VIBRATIONS

See (C15) in Appendix C. For this value ofis given byminimum

795

, the

(15)
from (C17) in Appendix C. By introducing the restriction, asolution quite different from (6) of the nonrestricted definitionin (4) is obtained. The second term of (15) corresponds tothe magnitude of the normalized cross-correlation function
and.
between
Let us denote , by whichtakes theminimum, that is(16)
during the period
For this , the phase change
determined and the average velocity in (3) is given by

is

(17)
where
denotes the phase value in the radian of.
by theBy multiplying the resultant velocity, the next object positionis estimated byperiodthe lag between the successively received complex signals is
uniquely determined.
The demodulated signal is A/D converted at a samplingof theinterval of . Nevertheless, the resultant estimate
next object position in (18) is represented not by a discretevalue, but by the continuous value in the above procedureis given from the
because the instantaneous velocityin (3). When the velocity is assumed
phase difference
to be 5 mm/s, for example, and the time intervalofs, the displacement
the transmission of the RF pulses isof the object ism during the time
, which is much less than the spatial resolution
interval
m of Fig. 1 in the depth directionis 1 s. Even if the samplingwhen the sampling period
period of 100 ns is employed, the resultant spatial resolution
m is still larger than the displacement ofm.
Thus, an accurate tracking of the object is realized by thismethod.
If the maximum value of the velocity is 0.1 m/s, theof the object during the time interval
displacement
ism, which is comparable with the spatial resolution of
m in the depth direction when the sampling
ns is employed. If the object position
period of, thisis on the boundary of the sampling interval of
displacement is not negligible and the possible values ofare, , andfor this case. For the case ofs, however, the spatial resolutionism and theabove-mentioned displacement of the object is neglected. For
this case, the proposed method does not differ from theconventional complex autocorrelation procedure except for therecursive adjustment of the sample window position based on
the displacement estimated by integration of the instantaneous
velocity.

(18)
From the resultant object position
and the velocity, the object position and the vibration velocity
on a large amplitude motion are simultaneously determined.
For the simple model described in (7) and (8), the minimum
is calculated as follows. By substituting,, andobtained in (B2)–(B6) of Appendix B intoof (C16) in Appendix
C, for the case(19)
and for the case
(20)
The characteristics of the minimumof (19) and (20)are shown in Fig. 3(b). As shown in this curve,
takes the minimum only at the true lag of. Thus, byintroducing the restriction into the cost function into (13),

D. Theoretical Comparison Between the Proposed
Method and the PLL-Based Method
In our method, by applying the constraint least-square
method to the demodulated reflected signal
and the preceding signal, the motion of the object isaccurately tracked from the phaseat the peak positionof the resultant complex correlation functions.
The proposed method differs from previous tracking methods in the following ways. In the standard cross-correlationbased method, as described in [11], it is well known thatthe velocity of the blood flow is obtained from the crosscorrelation if the blood vessel is at rest, that is, if the position of
the region of interest is unchanged. For the problem describedin this paper, however, there is a motion with large amplitude
and the next object position is estimated from the displacementof the object by multiplying the resultant velocity by theof the transmitted pulses. Thus, the spatial
time interval
resolution in the estimation of the object displacement isachieved in less than several micrometers in this paper, whichcorresponds to the cases where the reflected signal is A/D
converted with a sampling period of less than several ns (several micrometers/sound speed), which cannot be realized atpresent by the digitization of the data because a huge scaled

796IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

j^
j^ j
j^ j

j

Fig. 5. Theoretical comparison of error and the limitation vmax in theestimation of the velocity of an object by the previous PLL-based method
and the method proposed in this paper. The upper limit vmax of (21) and(22) is shown by the upper horizontal line. For the case where the absolutevalue of the velocity is less than the upper limit vmax , the velocity iscorrectly estimated by the proposed method as shown by the broken line.
In the PLL-based method, however, the object displacement during the timeinterval T is digitalized at a time interval =fCLK , that is, the spatial
resolution is =fCLK c0 = in the depth direction, as shown by the solidline. There is a quantization error which is not negligible.

1

Fig. 4. An illustration explaining theoretical limitations in the estimation of
the velocity of an object by 1) the previous PLL-based method and 2) the
method proposed in this paper: (a) for the signal at a time t, (b) an example
T , and (b0 ) another example of theof the succeeding signal at a time t
T.
succeeding signal at a time t+1+1

memory is necessary for continuous digitization for periods ofseveral heartbeats.
As shown in Fig. 4(2), based on the sampling theorem, theabsolute magnitude of the phase difference during the timemust be less than . Let us assume that theintervalismagnitude and the phase of the complex signal
represented by Fig. 4(2)(a) at a time . At the next time,of the case in Fig. 4(2)(b), the velocity is correctly determined,
while for the case Fig. 4(2)(b ) the phase differenceis equal to , which corresponds to the upper limit.
in
From (17), this condition determines the upper limit
the velocity estimation as follows:(21)
Fig. 5 explains the relationship between the true velocity andthe velocity estimates. The upper limitof (21) is shownby the upper horizontal line in Fig. 5. For the case where
the absolute value of the velocity is less than the upper limit, the velocity is correctly estimated by the proposed
method as shown by the broken line in the same figure.

1

2

2

1

On the other hand, let us consider the upper limit in
the previous PLL-based method. In the PLL-based method,
the zero-cross timing of the successively received RF signal
is searched for around the position which wasby adetermined for the previously received RF data
as shown indigital circuit with the clock frequency
Fig. 4(1).
isThus, if the object displacement during the interval
as shown in the relationship between
more thanFig. 4(1)(a) and (1)(b ), the sampling theorem does not hold.
Thus, the upper limit of this PLL-based method is described by(22)
which of course coincides with that of (21). In the PLL-based
method, however, the object displacement during the time
is digitalized at a time interval, that is,intervalin the depth direction,
the spatial resolution is
as shown by the solid line in Fig. 5. There is a quantization
error.
s,MHz,
For example whenm/s,s, andMHz, the upper limitin (21) and (22) is equal to 0.469 m/s. The quantization error in
the search for the zero-crossing point in the PLL-based method. This error
is equal to 0.188 m/s which is 40% that ofMHz.
is not negligible even for the caseAs described above for the PLL technique, the spatial
resolution of the tracking is limited to the clock frequency( several tens of megahertz) of the employed circuit. Thus,

KANAI et al.: HEART WALL VIBRATIONS

797

the PLL technique cannot reach the spatial resolution obtainedby the method proposed in this paper. The method proposed
herein makes it possible to continuously transcutaneously
measure the small vibration signals on a motion with a large
amplitude resulting from the heartbeat in a wide frequency
range.
III. A MEASUREMENT SYSTEM
In order to realize the procedure proposed in Section II-C,
a high-speed A/D converter with a large-scale memory and a,workstation are employed to analyze the complex signal,
resulting from the quadrature modulation of the signal
received by the sector-type ultrasonic transducer in the standard ultrasonic diagnostic equipment (Toshiba SSH-160A).
In the diagnostic equipment, the standard B-mode crosssectional image and M-mode image are displayed to identify
the measurement points on the heart wall. The frequencyof the transmitted ultrasound and repetition intervalofthe transmission of the ultrasonic pulses are employed as 3.75
MHz and 222 s, respectively, in the experiments detailed in
Sections IV, V, and VI.
, received by the ultrasonic transducer
The signal,is amplified and demodulated by the ultrasonic diagnostic, and imaginary
equipment. The resultant real signal,, are simultaneously A/D converted with a twosignal,. The
channel 12-b A/D converter at a sampling rate ofsampling rate in the system is variable from 1 to 10 MHz,and herein 1 MHz is employed in the following experiments.
In the A/D converter employed, the length of each signal
is limited to 1 046 000 points ( 1 megawords). However,by digitizing the real and imaginary signals only for thewhere the signal reflected at the object
period of length
is received by the ultrasonic transducer, the effective lengthof the digital signal in the A/D converter is increased fromone second to several seconds ( several beat periods). For
this intermittent digitization, the sampling clock of 1 MHz
by an external signal
is generalized only for the period
generator, which is completely synchronous with the masterclock of the ultrasonic diagnostic equipment. The demodulated
clock is also synchronous with the master clock.
Then, by applying the proposed method described in Section
II-C to the resultant digital signals, the velocity signal, denoted, is obtained at the workstation, which is connectedbywith the A/D converter via a GPIB interface. The electrocardiogram (ECG) and phonocardiogram (PCG) are also digitizedso that their timing is synchronous with that of the ultrasonicwave transmission from the transducer. The total time fromthe beginning of the A/D conversion to the display of the
waveform on the display is several minutes.
IV. VALIDATION EXPERIMENTS USING A WATER TANK
To confirm the principle of the proposed method, a smallof a rubber plate in a water tank asvibration signal
shown in Fig. 6 is estimated from the signal received by an
ultrasonic transducer. A motion with a large amplitude ( 7.5mm, 1 Hz), which corresponds to heartbeat, is generated by
the rotation of an eccentric cam. A vibration with a small

Fig. 6. An experimental system which simulates the vibration on the ventricle wall with a heartbeat of large amplitude using a vibrator and a motor withan eccentric cam. The motor with the eccentric cam generates motion having
a large amplitude of 7.5 mm and a frequency of 1 Hz which corresponds toa heartbeat. The vibrator imposes a small amplitude vibration ranging from30 to 1000 Hz with amplitude of 20 m peak to peak on this motion.

6

amplitude, which corresponds to the vibration on the wall of aventricle, is generated by a vibrator (Brüel & Kjær, 4810) on
the eccentric cam. In this experiment, the vibrator is excitedby band-limited noise ranging from 30 to 1000 Hz with an. The vibrator generates displacement of
amplitude of 5 V
the rubber plate which corresponds to the input voltage signal.
The amplitude of the corresponding displacement due to thevibrator is about 20 m peak to peak.
Fig. 7(a)–(c) shows the standard B-mode cross-sectional
image, the velocity, and M-mode images, respectively, which
are obtained by the standard ultrasonic diagnostic equipment.
As found in Fig. 7(b), the velocity of the rubber object cannot
be obtained for the whole period, since the measured point isfixed in space and the length of the sample volume is shorter
than the amplitude of the motion due to the eccentric cam. The
motion with large amplitude is shown in the M-mode image of
Fig. 7(c), however, the waveform of the small vibration due
to the generator cannot be recognized. One of the purposesof this paper is to measure the small vibration signal, whichshould be within the thickness of the curved line due to the
motion of Fig. 7(c) in the M-mode image.
By applying the proposed method in Section II-C to thedigital signals in the measurement system described in Section, is obtained. The
III, the velocity signal, denoted by
input voltage signalof the vibrator and the output signalof a laser Doppler velocimeter (Ono Sokki LV1300,
Hakusan 1-16-1, Midori-ku, Yokohama, Japan) are also A/Dconverted simultaneously.
Fig. 8 shows the results of these experiments. Fig. 8(b)of the vibrator. Fig. 8(a) shows
shows the input signals
obtained by thethe estimates of the object motion
proposed method in (18), which is superimposed on the Mmode image generated from the magnitude of the demodulated. In these experiments, since the length
received signalof the transmitted RF signal is about 3 s, the range

798IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

Fig. 7. Results obtained by standard ultrasonic equipment in the experiments on the detection of the vibration signals on the rubber plate of the vibrator
on the eccentric cam in the water tank of Fig. 6. (a) B-mode image showing the cross-sectional image around measured point A: The scale of the verticaland that of the horizontal axes are 1 cm, and (b) the velocity measured at a fixed position of the sample volume, which is indicated at point A in (a).
The scales of the vertical axis and the horizontal axis are 0.05 m/s and 100 ms, respectively, and (c) the M-mode image showing the slow motion on
the axis across point A. The scale of the vertical axis is 1 cm in depth.

employed in the calculating correlation function in (13) is 1
point in length. The tracking resultof point exactly
coincides with the M-mode image as shown in Fig. 8(a).
Fig. 8(d) shows the vibration velocity signal
estimated by the proposed method in (17). The amplitude of
the signal
in Fig. 8(d) is similar to that of the signals
obtained by laser Doppler velocimeter in Fig. 8(c). A
high-pass filter with a cutoff frequency of 0.5 Hz is employedin the laser Doppler velocimeter. Due to this filter, it isexperimentally found that there is a time delay of about 850
ms only for the low-frequency component of 1 Hz generated
by the rotation of the eccentric cam. By considering thetime delay, these waveforms ofandalmostcompletely coincide in the low-frequency components.
To compare the high-frequency components of these signals,
Fig. 9(a)–(c) shows the waveforms of the input voltage signal

KANAI et al.: HEART WALL VIBRATIONS

799

Fig. 8. Estimation results using the experimental system shown in Fig. 6:(a) the object motion x t estimated by the proposed method, superimposedon the M-mode image, obtained by the magnitude of the received signal,(b) input voltage signal Vinput t of the vibrator by which a vibration of
band-limited noise ranging from 30 to 1000 Hz with an amplitude of 20 mpeak to peak is generated, (c) the velocity signal vlaser t obtained by laser
Doppler velocimeter, and (d) the vibration velocity signal vus t estimatedby the proposed method.

^( )

()

()

^ ()

Fig. 9. To compare the high-frequency components of the waveforms in
Fig. 8, (a)–(c) show the waveforms of the input voltage signal Vin t of
the vibrator, the output signal vlaser t of a laser Doppler velocimeter,
and the vibration velocity signal vus t estimated by the proposed method,respectively, during the period from 6.27 to 6.33 s in Fig. 8.

()

()
^ ()

of the vibrator, the output signalof a laser
Doppler velocimeter, and the vibration velocity signal
estimated by the proposed method, respectively, during theperiod from 6.27 to 6.33 s in Fig. 8. Though there is a slightdifference in amplitude of these signals, their waveforms are
similar. The result in Fig. 9(b) obtained by the laser Doppler
is correct as to amplitude. In the proposed procedure using
the ultrasonic Doppler, however, the instantaneous velocity
is obtained from (3) by assuming that there is only one
frequency component with a center frequency component ofin the employed ultrasonic wave. However, thepulsive ultrasonic wave has a frequency band of about 1
MHz around the center frequency. At the same time, there islarger attenuation for the higher frequency components of theultrasonic wave. Thus, the velocity in Fig. 9(c) obtained by theultrasonic Doppler-based method is slightly underestimated.
To quantitatively evaluate the correlation at each freamong the estimated signal
in Fig. 8(d),quencyin Fig. 8(c), andin Fig. 8(b), the following
andsquared magnitudes of coherence functions
are evaluated:
(23)
(24)
where

,

, and
are the spectra of,, and, respectively. The term
denotes theaverage operation. There are 27 averaging operations in thisexperiment. The squared magnitude of coherence function
quantitatively evaluates the correlation between the two signals
for each frequency component [33]. If the squared magnitude, between two signals is equalof coherence function,
to one, there is complete correlation between these signals forthe frequency component . That is, one signal is completely$ ()Fig. 10. (a) The squared magnitude of coherence function j in us f j2
between Vin t and vus t , (b) the magnitude of the transfer function
Hin us f from the input signal Vin t to vus t estimated by the proposed
method in the frequency range up to near the Nyquist frequency = T , and(c) the phase characteristics of the transfer function Hin us f .

! ()

()

^ ()

() ^ ()

11! ()described by the other signal and the linear relationshipbetween them.
fromAt the same time, the transfer functiontoand the transfer function
from

800

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996$

Fig. 11. (a) The squared magnitude of coherence function j laser us (f )j2
between vlaser (t) and v^us (t), (b) the magnitude of the transfer function
Hlaser us (f ) from vlaser (t) to v^us (t) estimated by the proposed method,and (c) the phase characteristics of the transfer function Hlaser us (f ).

!

!

to

$
$

6

are obtained as follows [33]:$Fig. 12. The average values 0in us and 0laser us of the squared magnitude of coherence functions j in us (f )j2 and j laser us (f )j2 , respectively, for various values of 1x of the range R = [x(t) 1x; x(t) + 1x]
employed in (13). (a) The average values from the frequency ranges from 1
Hz to 1 kHz and (b) the average values from the frequency ranges from 900
Hz to 1 kHz. The optimum value of 1x is 1 point in length, that is, therange R corresponds to about 0.75 mm around the previously determined
object position x(t).

$

0

6

(25)
(26)
Fig. 10(a)–(c), respectively, shows the resultant squared
magnitude of coherence function
between, and, and the magnitude and the phase of thefrom the input signal
totransfer functionin the frequency range up to near the Nyquist frequency. Fig. 11(a)–(c), respectively, shows the resultant
squared magnitude of coherence function
betweenand, and the magnitude and thefromphase of the transfer function
to.
From Figs. 10(a) and 11(a), the squared magnitude ofcoherence function is almost one from 30 Hz to 1 kHz. Thereandare, therefore, linear relationships between
and between
and. The input signal
does not include low-frequency components less than 30 Hz,, the low-frequency components less than 100
while forHz resulting from rotation of the eccentric cam are included.
in Fig. 10(a) is almost zero in the lowThus,
frequency range less than 100 Hz, in which low-frequency
movements due to the eccentric cam are not negligible. However, since there are low-frequency components with large

amplitude due to the eccentric cam in bothand,in Fig. 11(a) is almost one even in
the low-frequency range.
From Fig. 10(c), the phase characteristics are described byan almost straight line in the frequency range from 100 Hz
to 1 kHz. Though the phase characteristics
in the low-frequency components of Fig. 11(c) are affectedby the high-pass filter of 0.5 Hz employed in the laser
Doppler velocimeter, the magnitude of the transfer functionof Fig. 11(b) is almost flat in the frequency
range from 1 Hz to 1 kHz. There are differences from theoptimum characteristics (flat magnitude and linear phase)of the transfer function. However, these differences can becorrected because the squared magnitude of coherence isalmost one in the frequency range less than 1 kHz. From theseexperimental results, it is confirmed that the small vibrationon the large motion is successfully detected by the proposed
method in the frequency range from 1 to 1000 Hz.
Fig. 12 shows the average values, respectively denoted byand, of the squared magnitude of coherencefunctionsandfor various valuesofof the range
employedinof (13). Fig. 12(a) and (b) shows the averagevalues for the frequency ranges from 1 Hz to 1 kHz and from

KANAI et al.: HEART WALL VIBRATIONS

801

Fig. 13. A standard B-mode image showing the cross-sectional area around the two measured points, (A) and (B ), in an in vivo experiment for the detectionof vibrations on the interventricular septum of the left ventricle (LV) of a 23-year-old male patient with acute lymphoblastic leukemia who has been treatedwith antracenadiones. Points (A) and (B ) are on the RV and LV surfaces of the interventricular septum, respectively. The ultrasonic beam passing throughthe two points is almost perpendicular to the septum during the measurements.

900 Hz to 1 kHz, respectively. For the frequency range from900 Hz to 1 kHz in Fig. 12(b), the results vary remarkably.
, the range
where the constraint
For the case ofcross-correlation function is evaluated has only one point andthe velocity estimates are easily affected by additive noise.
increases, the effect of noise is decreased. At
As range
the same time, the spatial resolution is also decreased. Fromis 1
Fig. 12, it is found that the optimum value ofcorresponds to 0.75 mm
point in length, that is, range. Thisaround the previously determined object position
value of1 is employed in the following in vivo
experiments.
V. IN VIVO DETECTION OF HEART WALL VIBRATIONSIn order to noninvasively diagnose myocardial damage byadriamycin injection based on the acoustic characteristics ofthe myocardium, let us try to measure the small vibration
signals on the heart wall from the precordium and apply spectrum analysis to the resultant vibration signals. It is, however,difficult to measure a small vibration signal with an amplitudeless than 100 m with the previously developed methods
because such a small vibration is superimposed on the large
motion resulting from a heartbeat with an amplitude of about
10 mm. For this problem, we apply the proposed method
for noninvasive measurement of small vibration signals on theinterventricular septum as follows.
First, the proposed method is applied to the detectionof a vibration signal on the interventricular septum of a23-year-old male patient with acute lymphoblastic leukemia
who has been treated with antracenadiones. Fig. 13 shows
the B-mode image around two measured points, ( ) and( ). This image is obtained by standard ultrasonic diagnostic
equipment. Points ( ) and ( ) are on the surface of the rightventricular (RV) side and the left ventricular (LV) side of
the interventricular septum, respectively. The ultrasonic beampassing through points ( ) and ( ) is almost perpendicular to
the interventricular septum as shown in Fig. 13.
Fig. 14(b) and (c) shows the ECG and the heart sound(PCG), respectively. The tracking results
and,estimated by (18) of the movement of points ( ) and ( )are superimposed on the M-mode image in Fig. 14(a), which
was obtained from the magnitude of the demodulated received
signals. Fig. 14(d) and (e) shows the estimates of the vibration
velocity signals
andof the tracked points
and, respectively. The resultant vibration signals are
sufficiently reproducible for the six heartbeat periods as shownin Fig. 14(d) and (e).
For the systolic phase, where the ventricular contraction
begins near the peak of the R-wave of the ECG, the velocity
estimates
andof the septum have positive values
after negative steep peaks as shown in Fig. 14(d) and (e); thatis, the septum moves in the opposite direction of the ultrasonictransducer on the chest wall. The amplitude ofon theLV side is a little greater than that ofon the RV side,
especially in the systolic phase. This is due to the thickness
change of the interventricular septum.

802IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

Fig. 15. The waveforms of ECG, PCG, v^A (t), and v^B (t) of Fig. 14 froma P-wave to the next P-wave in the succeeding heartbeat are separated. The
results for the first five heartbeats in Fig. 14 are overlaid and shown: (a) the
ECG, (b) the PCG, (c) and (d) the estimates of the vibration velocity signalsv^A (t) and v^B (t) of points A and B . The reproducibility of v^A (t) and v^B (t)is qualitatively confirmed at least for the low-frequency components.

Appendix D:
Fig. 14. In vivo experimental results of the vibration at points (A) and (B )on the RV and the LV surfaces of the interventricular septum, respectively,
in Fig. 13: (a) the tracking results x^A (t) and x^B (t) in (18) of the points A
and B . These results are superimposed on the M-mode image, (b) the ECG,(c) the PCG, (d), (e) the estimates of the vibration velocity signals v^A (t)^and v^B (t) of points A and B, and (f) the change in thickness, hAB (t), of^
the interventricular septum given by (27). The waveform of hAB (t) is also
reproducible for the six heartbeat periods.

Fig. 14(f) shows the thickness change,, of theinterventricular septum from the thickness at the timing ofis given by
the R-wave of ECG.

(27)

The waveform ofis also reproducible for the six
heartbeat periods. For the systolic phase, the interventricular
septum becomes about 1 mm thicker than that of the diastolic
phase.
, andof Fig. 14
The waveforms of ECG, PCG,
from a P-wave to the next P-wave in the succeeding heartbeat
are separated. The results for the first five heartbeats in Fig. 15
are overlaid and shown in Fig. 14. From Fig. 15(c) and (d), theandis qualitatively confirmed
reproducibility ofat least for the low-frequency components less than 20 Hzwith large amplitude resulting from the heartbeat.
To quantitatively evaluate the reproducibility of the simsignals,
forilarity of theeach frequency component of , the following reproducibilis introduced by (D5) as described in
ity function(28)

is the Fourier spectrum of. Similar to thewheresquared magnitude of the coherence function,also
is equal to one,has a value from zero to one. If
waveforms
completely coincide for the frequency
is zero, there is no correlation
component of . If
amongfor the frequency components of .
The solid lines in Fig. 16(a) and (b) show the reproducibility, evaluated for the components around thefunction,
andofsecond heart sounds of the vibration signals
the first five heartbeats in Fig. 14(d) and (e), respectively. Eachspectrum is obtained by applying the fast Fourier transformandfor the period of 100 ms around the(FFT) tosecond heart sound of each heartbeat. For these results, there
is a high correlation between these small vibration signals forfive heartbeats up to at least 100 Hz. Thus, the reproducibilityof the measured vibrations is quantitatively confirmed.
The solid lines in Fig. 16(c) and (d) show the average
andof the first five
power spectra of the signals
heartbeats in Fig. 14(d) and (e), respectively. The analyzed
periods are same as those in Fig. 16(a) and (b). The minimum
and maximum values of the power spectra of the vibration
signals for five heartbeats are shown by the vertical bars foreach frequency. As shown in these figures, there are large
components in the low-frequency range less than 25 Hz, whichcorrespond to the motion due to the heartbeats. However,there is enough reproducibility even for the high-frequency
components from 25 to 90 Hz. The frequency component at
90 Hz is less than the maximum component at 10 Hz by

KANAI et al.: HEART WALL VIBRATIONS

Fig. 16.
j803(a) and (b) The solid lines show the reproducibility function,

2^A (t) and v^B (t) for the first five0 (f )j , of the small vibration signals v

heartbeats in Fig. 14(d) and (e), respectively. Each spectrum is obtainedby applying FFT to v^A (t) and v^B (t) for the period of 100 ms aroundthe second heart sound for each heartbeat, where the Hanning window isemployed. There is a high correlation between these small vibration signals
for five heartbeats up to at least 100 Hz. The broken lines show the squared
magnitude of the standard coherence function between v^A (t) and PCG andbetween v^B (t) and PCG, respectively. The small vibration detected on theinterventricular septum at the timing around the second heart sound has a highcorrelation with the second heart sound. (c) and (d) The average power spectraof v^A (t) and v^B (t), respectively. The analyzed periods are the same as thosein (a) and (b). The vertical bar for each frequency shows the minimum andmaximum values of the power spectra of the vibration signals.

6Fig. 17.

The squared magnitude of the standard coherence function,j AB (f )j2 , between v^A (t) of Fig. 14(d) and v^B (t) of Fig. 14(e) at theperiod 6100 ms around the second heart sound. There is a correlation

between these small vibration signals of both sides of the interventricular
septum up to at least 100 Hz.

about 25 dB; that is, the amplitude of the components at
90 Hz is 1/20 that of the components at 10 Hz. For thesesmall vibration components, which are within the thickness

fgFig. 18. (a) and (b) PCG and ECG measured when vB (t) of (c) are
obtained. (c) and (d) The vibrations vB (t) at the LV side of the interventricular septum of a presumedly healthy 22-year-old male volunteer,where the measured point B is selected in the cross-sectional image along
the longitudinal axis and the short axis, respectively, of the heart. The
small vibration signals of six heartbeats are overlaid in each figure and the
waveforms are shown for the period around the second heart sound. (e) and (f)Average power spectra of the signals in (c) and (d), respectively. The average
power spectra of the vibration signals for six heartbeats are shown by thesolid lines. The minimum and maximum values of the vibration signals for
six heartbeats are shown by the vertical bars for each frequency component
in the power spectra. The waveforms in (c) are similar to those in (d). The
power spectrum in (e) is also similar to that in (f).

fgof the overlaid curves in Fig. 15(c) and (d), the waveform isaccurately measured by the proposed method.
The broken lines in Fig. 16(a) and (b) show the squared
magnitude of the standard coherence function between
and PCG and between
and PCG, respectively. It is
known that the second heart sound is generated due to theclosure of the aortic valve [34]. The small vibration on theinterventricular septum measured here has a high correlation
with the second heart sound, which is a clue for identifying
the origin of the heart wall vibration detected here.
Fig. 17 shows the squared magnitude of the standard coherence function,, between
andat theperiod 100 ms around the second heart sound. There is a804IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

Fig. 19. The waveforms of the vibration signal on the LV side of the interventricular septum of three presumedly healthy young male volunteers and theiraverage power spectra. (a) A presumedly healthy 22-year-old male volunteer; eight heartbeats overlaid, (b) a presumedly healthy 26-year-old male volunteer;
five heartbeats overlaid, and (c) a presumedly healthy 22-year-old male volunteer; six heartbeats overlaid. The waveforms of the several heartbeats are overlaid
and are shown at the moment of the second heart sound emission. For each subject, the reproducibility is quantitatively confirmed up to 100 Hz. For the power
spectra, the employed liens are shown by a format similar to that in Fig. 18(e) and (f). The waveforms and power spectra of these volunteers are similar.

correlation up to at least 100 Hz between these small vibration
signals on both sides of the interventricular septum.
As shown in Fig. 13, measured points ( ) and ( ) on
the interventricular septum are selected in the cross-sectional
image along the longitudinal axis of the heart. The direction of
the ultrasonic beam passing through the measurement points
is selected to be almost perpendicular to the interventricular
septum during the measurements. There may be, however,complex movement in the heart wall during the period of
one heartbeat. Although the ultrasonic beam is perpendicular
to the septum during the measurement in the cross-sectional
image along the longitudinal axis, it is not guaranteed that theultrasonic beam is perpendicular to the wall surface because
it is probable that the ultrasonic beam tilts in the crosssectional image along the short axis. Thus, the small vibrationis also measured for the point selected in the cross-sectional
image along the short axis of the heart and is comparedwith that for the longitudinal axis as follows. Fig. 18(c) and(d) shows the small vibrations
at the LV side of
the interventricular septum of a presumedly healthy 22-yearold male volunteer, where the measured point B is selectedin the cross-sectional image along the longitudinal axis andthe short axis, respectively. The small vibration signals of
six heartbeats are overlaid in each figure and the waveforms
are shown for the period around the second heart sound in
Fig. 18(a). Fig. 18(e) and (f) shows the average power spectraof the signals in Fig. 18(c) and (d), respectively. The average
power spectra of the vibration signals for six heartbeats areshown by the solid lines. The minimum and maximum valuesof the vibration signals for six heartbeats are shown by thevertical bars for each frequency component in the power
spectra. For the frequency components less than 80 Hz, thevertical bar is very short, that is, sufficient reproducibility is
obtained. The waveforms in Fig. 18(c) are similar to those in
Fig. 18(d). At the same time, the power spectrum in Fig. 18(e)is also similar to that in Fig. 18(f). From these results, it isfound that the ultrasonic beam is perpendicular to the wall
plane of the interventricular septum during the measurements.
Thus, the small vibration signal of the interventricular septumis measured with sufficient reproducibility by selecting theultrasonic beam passing through the measured points so that it
is perpendicular to the heart wall in the cross-sectional images
along the conventional longitudinal axis or the conventional
short axis of the heart.
VI. COMPARISON OF HEART WALL VIBRATIONSBETWEEN NORMAL SUBJECTS AND PATIENTS
Fig. 19 shows the waveforms of the small vibration signals on the LV side of the interventricular septum of threepresumedly healthy young male volunteers and their average power spectra. The waveforms of several heartbeats areoverlaid and shown at the moment of the emission of thesecond heart sound. The horizontal center in Fig. 19(a-1)–(c-1)
corresponds to the moment of the emission of the second heartsound. For each subject, the reproducibility is quantitatively
confirmed up to at least 100 Hz using the reproducibilityof (28). For the power spectra, the employedfunctionlines are the same format as those in Fig. 18(e) and (f). By
comparing Fig. 19(a-1)–(c-1) and (a-2)–(c-2), it can be seenthat the waveforms and power spectra of these three volunteers
are similar.
For three patients with acute lymphoblastic leukemia whohave been treated with antracenadiones (mitoxantrone), thewaveforms of the vibration signals on the LV side of theKANAI et al.: HEART WALL VIBRATIONS

805

Fig. 20. The waveforms of the small vibration signals on the LV side of the interventricular septum for three patients with acute lymphoblastic leukemia
who have been treated with antracenadiones (mitoxantrone), and their average power spectra. (a) A 32-year-old male patient; eight heartbeats overlaid; themeasurement was performed two months before his death. (b) The same patient as (a); eight heartbeats overlaid; the measurement was performed three months
before the measurement of (a). (c) A 23-year-old male patient; five heartbeats overlaid; this subject is also represented in Figs. 13–17. (d) A 25-year-old
female patient; six heartbeats overlaid. For the waveforms of these patients, high-frequency components higher than 50 Hz are decreased.

interventricular septum and their average power spectra areshown in Fig. 20. The subject in Fig. 20(c) is the same asthat represented in Figs. 13–17. The subjects in Fig. 20(c) and(d) are in complete remission. The subject in Fig. 20(b) wasin the hospital. Three months after the measurement of the
waveforms in Fig. 20(b), the waveforms in Fig. 20(a) were
measured for the same patient. Then about two months after
the last measurement, the subject died. For the four waveformsof these patients, the frequency components higher than 50
Hz are less as compared with the power spectra in Fig. 19.
Especially for the waveforms in Fig. 20(a-1), only the lowfrequency components less than 50 Hz are included. For
the subject in Fig. 20(d-1), the components lower than 30
Hz are decreased. For the subject in Fig. 20(c-1), the poweris still large for the frequency components up to 80 Hz.
However, the waveforms during the isovolmetric relaxation
period, which begins just after the moment of the second
heart sound emission, are different from those for the healthy
volunteers in Fig. 19.
To quantitatively evaluate the decrease of the vibration
power of the patients from those of the healthy volunteers,
three kinds of average power,dB from dc to 30 Hz,
dB from 30 to 60 Hz, anddB from 60 to 100 Hz, arecalculated from the average power spectra in Figs. 19 and 20.
The distribution of the resultant values of the average power,, , and , is indicated by circles for the three healthy
subjects and by x’s for the patients in Fig. 21. As shown in

Fig. 21. Distribution of the three values of the average power, P1 dB fromdc to 30 Hz, P2 dB from 30 to 60 Hz, and P3 dB from 60 to 100 Hz foreach patient, denoted by x’s, of Fig. 20 and each healthy volunteer, denotedby circles, of Fig. 19. The power of the heart wall vibrations of patients isdecreased by at least several dB in comparison with those of the healthy
subjects.

Fig. 21, the power of the vibration of all patients is decreased
by at least several dB as compared with that of the healthy
subjects.
Myocardial physical properties have been reported tochange sensitively under various pathological conditions. In
a study using an isolated canine left ventricle with regional
ischemia, the instantaneous impedance frequency curve [35]changed from a single-peak configuration to a double-peak
configuration after the cessation of coronary flow and returnedto the single-peak configuration with the recommencement of806IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

coronary flow [36]. The difference between the power spectrain Fig. 19 and those in Fig. 20 corresponds to the changes in
the acoustical characteristics of the local myocardium dueto the injection of mitoxantrone. Mitoxantrone has beenreported to have an adverse cardiotoxic effect, which resultsin irreversible diffuse myocardial fibrosis at cumulative or
high doses.
Taking these experimental results into account, the method
proposed in this paper provides a new approach to the detection of such myocardial physical and/or histological heterogeneity of a ventricle by analyzing human heart vibrations
measured noninvasively from the chest surface.
VII. CONCLUSIONSIn this paper we have proposed a new method for accurate
measurement of small vibrations on a large motion using ultrasound. Using a measurement system to realize this method,
small vibrations were accurately measured in the frequency
range up to 1 kHz in experiments using a water tank. By
applying the method and the measurement system to the measurement of vibration signals of the interventricular septum inin vivo experiments, sufficient reproducibility was confirmed
for the frequency components up to 100 Hz. With the systemdevised herein, it is possible to detect small vibrations up to 1
kHz. However, the amplitude of frequency components higher
than 100 Hz in the in vivo experiment was too small to bemeasured accurately. The obtained change in thickness of the
septum will directly relate to the regional contractility of the
myocardium and it is diagnostic for ischemia and infarction.
Discussion of the data on thickening and their implications
will be presented in a separate paper in the near future.
Further investigation of this proposed method, including
its clinical application to the noninvasive local diagnosis ofcoronary artery disease, drug-induced myocardial disease, andarteriosclerosis, is being conducted. It is necessary to identify
the origin of and the process by which the heart wall vibrations
detected in this paper were generated. In this paper, measured
vibration signals are analyzed in the frequency domain only forthe period around the second heart sound emission. It is also
necessary to analyze the vibration signals for other periods.
APPENDIX ANONCONSTRAINT LEAST-SQUARE APPROACH
This appendix describes the minimization of the normalizedof (4). Equation (4) ismean squared difference
simply rewritten as follows:

and
are a
where denotes the complex conjugate,
real function of and a real constant, respectively, defined by
(A2)
(A3)
and

is a complex function ofdefined by(A4)

Describing
andby the sum of their real andand, respectively,
imaginary parts,
the normalized mean squared differencein (A1),, is described bymultiplied by(A5)
By taking partial derivatives ofandand setting them to zerowith respect to

(A6)

(A7)
From (A6) and (A7)
(A8)
Thus
(A9)
By substitutingofof (A9) into (A1), the minimum value
is given by(A10)
Using the original signals
minimum value
as follows:

and
of

, the
in (A1) is obtained

(A1)

(A11)
This equation is inferred by (6) of Section II-B.

KANAI et al.: HEART WALL VIBRATIONS807

APPENDIX BDERIVATION OF THE VALUES OF THE CONSTANTS

For the case of
,

, AND

This appendix theoretically derives the values ofofof (A3) and (C3), andof (A4) and(A2) and (C2),(C4) for the simple model of (7) and (8) in Section II-B. The
defined in (A2) and (C2) is obtained byfunction(B6)
while for the case(B1)
where the range

is defined by(B7). For the case
APPENDIX C
A CONSTRAINT LEAST-SQUARE APPROACH

This appendix describes the minimization of the normalizedof (13) under the restriction
mean squared difference
that there is only a phase change between the signals
and. Equation (13) is simply rewritten as follows:(B2)
while for the case(B3)
For the constant(C1)

defined in (A3) and (C3)
where
and
are a real function of
constant, respectively, defined by

and a real
(C2)
(C3)

and
(B4)
The complex functionbyis a complex function ofdefined by(C4)

defined in (A4) and (C4) is given
Describing
and by the sum of their real and imaginary
and, respectively, the normalizedparts,in (C1), multiplied bymean squared difference, is described by

(B5)

(C5)

808

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996

By taking partial derivatives of, and and setting them to zerowith respect to,ofinto (C1), the minimum value
is given byin (C1)

(C6)

(C7)

(C8)
From (C6) and (C7)(C16)

Using the original signalsofminimum value

and

, the
in (15) is given by

(C9)
Thus
(C10)

By substituting (C10) into (C8)
(C11)

(C17)
This equation is inferred by (15) of Section II-D.
There is another procedure for this minimization of theof (C1) undernormalized mean squared difference
as follows. By setting to,
the restriction of(C1) is described by the nonrestricted form as

becauseandare real constants. Thus, the following
quadratic equation is obtained:
(C12)

(C18)

The solutions are given by(C13)
By substituting these solutions ofvalues of are obtained by

By taking partial derivatives ofsetting it to zerowith respect toandinto (C10), the optimum(C19)
By multiplying (C19) by
(C14)

(C20)
Thus, the optimum value ofis given by

(C15)
where

denotes the phase of the complex function. From (C4) and (C15), (14) is obtained.
of (C1) is
For , the mean squared difference
maximized, while for, the minimization ofisof (C15) and the relation of (C8)achieved. By substituting

(C21)

(D1)

KANAI et al.: HEART WALL VIBRATIONS

809

These two solutions are identical to those in (C14) and (C15).
The form of this procedure is direct. In this paper, however,the constraint is explicitly described by introducing restrictionof the second term into the right-hand side of (C1).
APPENDIX D
QUANTITATIVE EVALUATION OF REPRODUCIBILITY
signals
To quantitatively evaluate the similarity of the,() for each frequency component ,of the spectrum
the normalized mean squared differenceoffrom the spectrumof a waveform
is defined in (D1), as shown at the bottom of the preceding
denotes the number of heartbeats. By taking
page, wherewith respect toa partial derivative of(D2)
Thus, the optimum spectrum
is given by, which minimizes,(D3)
which coincides with the average of
By substituting (D2) and the optimumof (D1)

as a natural result.
of (D3) into

(D4)

Let us define the reproducibility function,
second term of the equation by, from the

(D5)

Since
has values from zero to one, the value ofis also from zero to one. If
is equal to one,
waveforms
completely coincide for the frequency
is zero, there is no correlation
component of . If
for the frequency components of .
between

ACKNOWLEDGMENT
The authors would like to thank Prof. Emeritus M. Tanaka of
Tohoku University, Prof. F. Dunn of the Bioacoustics ResearchLaboratory of the University of Illinois, and Prof. K. Shirado,
Assoc. Prof. Yoshiko Saito, and Dr. E. Kamata of Tohoku
University, School of Medicine, for their suggestions. The
authors also wish to acknowledge the contributions of M.
Takano in their laboratory.
REFERENCES[1] H. Kanai, N. Chubachi, K. Kido, Y. Koiwa, T. Takagi, J. Kikuchi, and T.
Takishima, “A new approach to time-dependent AR modeling of signals
and its application to analysis of the fourth heart sound,” IEEE Trans.
Signal Processing, vol. 40, pp. 1198–1205, May 1992.
[2] S. Satomura, “Ultrasonic Doppler method for the inspection of cardiac
function,” J. Amer. Soc. Acoust., vol. 29, no. 11, pp. 1181–1185, Nov.
1957.
[3] D. W. Baker, “Pulsed ultrasonic Doppler blood-flow sensing,” IEEE
Trans. Sonics Ultrason., vol. SU-17, pp. 170–185, July 1970.
[4] F. D. McLeod and M. Anliker, “A multiple-gate pulsed directional
Doppler flowmeter,” in Proc. IEEE Ultrason. Symp., Miami, FL, Dec.
1971.
[5] S. L. Johnson, D. W. Baker, R. A. Lute, and H. T. Dodge, “Doppler
echocardiography,” Circulation, vol. XLVIII, pp. 810–822, Oct. 1973.
[6] F. E. Barber, D. W. Baker, A. W. C. Nation, D. E. Strandness, Jr.,and J. M. Reid, “Ultrasonic duplex echo-Doppler scanner,” IEEE Trans.
Biomed. Eng., vol. BME-21, pp. 109–113, Mar. 1974.
[7] C. J. Hartley, H. G. Hanley, R. M. Lewis, and J. S. Cole, “Synchronizedpulsed Doppler blood flow and ultrasonic dimension measurement inconscious dogs,” Ultrasound Med. Biol., vol. 4, pp. 99–110, 1978.
[8] M. Brandestini, “Topoflow—A digital full range Doppler velocity meter,” IEEE Trans. Sonics Ultrason., vol. SU-25, pp. 287–293, Sept.
1978.
[9] E. Wildi, J. W. Knutti, H. V. Allen, and J. D. Meindl, “Dynamics andlimitations of blood/muscle interface detection using Doppler power
returns,” IEEE Trans. Biomed. Eng., vol. BME-27, pp. 565–573, Oct.
1980.
[10] W. D. Barber, J. W. Eberhard, and S. G. Karr, “A new time domain
technique for velocity measurements using Doppler ultrasound,” IEEE
Trans. Biomed. Eng., vol. BME-32, pp. 213–229, Mar. 1985.
[11] C. Kasai, K. Namekawa, A. Koyama, and R. Omoto, “Real-time twodimensional blood flow imaging using an autocorrelation technique,”IEEE Trans. Sonics Ultrason., vol. SU-32, pp. 458–463, May 1985.
[12] R. M. Olson and D. K. Shelton, “A nondestructive technique to measurewall displacement in the thoracic aorta,” J. Appl. Phys., vol. 32, no. 1,
pp. 147–151, Jan. 1972.
[13] R. M. Olson and J. P. Cooke, “A nondestructive ultrasonic technique
to measure diameter and blood flow in arteries,” IEEE Trans. Biomed.
Eng., vol. BME-21, pp. 168–171, 1974.
[14] J. O. Arndt, “The diameter of the intact carotid artery in man and itschange with pulse pressure,” Pflüegers Arch., vol. 301, pp. 230–240,
1968.
[15] D. J. Mozersky, D. S. Sumner, D. E. Hokanson, and D. E. Strandness,
Jr., “Transcutaneous measurement of the elastic properties of the humanfemoral artery,” Circulation, vol. XLVI, pp. 948–955, Nov. 1972.
[16] A. P. G. Hoeks, C. J. Ruissen, P. Hick, and R. S. Reneman, “Transcutaneous detection of relative changes in artery diameter,” UltrasoundMed. Biol., vol. 11, no. 1, pp. 51–59, 1985.
[17] W. T. Kemmerer, R. W. Ware, H. F. Stegall, J. L. Morgan, and R.
Kirby, “Blood pressure measurement by Doppler ultrasonic detectionof arterial wall motion,” Surgery, Gynecology, Obstetrics, vol. 131, pp.
1141–1147, Dec. 1970.
[18] K. Lindström, K. Marsal, G. Gennser, L. Bengtsson, M. Benthin, and
P. Dahl, “Device for measurement of fetal breathing movements, 1.
The TD-recorder. A new system for recording the distance between two
echogenerating structures as a function of time,” Ultrasound Med. Biol.,vol. 3, pp. 143–151, 1977.
[19] D. E. Hokanson, D. E. Strandness, Jr., and C. W. Miller, “An echotracking system for recording arterial wall motion,” IEEE Trans. Sonics
Ultrason., vol. SU-17, pp. 130–132, July 1970.
[20] D. E. Hokanson, D. J. Monzersky, S. D. Sumner, and D. E. Strandness,
Jr., “A phase-locked echo tracking system for recording arterial diameterchanges in vivo,” J. Appl. Phys., vol. 32, no. 5, pp. 728–733, May 1972.

810

IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRICS, AND FREQUENCY CONTROL, VOL. 43, NO. 5, SEPTEMBER 1996[21] K. Nakayama and S. Sato, “Ultrasonic measurement of arterial wall
movement utilizing phase-tracking system,” in Dig. 10th Int. Conf. Med.
Biol. Eng., Dresden, 1973, p. 318.
[22] D. N. White and R. J. Stevenson, “Transient variations in the systolicpulsations in amplitude of intracranial echoes, their artificial origin,”
Neurology, vol. 26, pp. 683–689, 1976.
[23] C. F. Olsen, “Doppler ultrasound: A technique for obtaining arterial
wall motion parameters,” IEEE Trans. Sonics Ultrason., vol. SU-24, pp.
354–358, June 1977.
[24] L. W. Korba, R. S. C. Cobbold, and A. J. Cousin, “An ultrasonic imaging
and differential measurement system for the study of fetal respiratory
movements,” Ultrasound Med. Biol., vol. 5, pp. 139–149, 1979.
[25] I. Rapoport and A. J. Cousin, “New phase-lock tracking instrument
for foetal breathing monitoring,” Med. Biol. Eng. Comput., vol. 20, pp.
1–6, Jan. 1982.
[26] D. H. Groves, T. Powalowski, and D. N. White, “A digital technique
for tracking moving interfaces,” Ultrasound Med. Biol., vol. 8, no. 2,
pp. 185–190, 1982.
[27] C. J. Hartley, H. Litowitz, R. S. Rabinovitz, W. X. Zhu, J. E. Chelley,
L. H. Michael, and R. Bolli, “An ultrasonic method for measuring tissue
displacement: Technical details and validation for measuring myocardial
thickening,” IEEE Trans. Biomed. Eng., vol. 38, pp. 735–747, Aug.
1991.
[28] L. S. Wilson and D. E. Robinson, “Ultrasonic measurement of small
displacements and deformations of tissue,” Ultrason. Imag., vol. 4, pp.
71–82, 1982.
[29] C. J. Hartley, L. A. Latson, L. H. Michael, C. L. Seidel, R. M. Lewis,
and M. L. Entman, “Doppler measurement of myocardial thickening
with single epicardial transducer,” Amer. J. Physiol., vol. 245, pp.
H1066–1072, 1983.
[30] A. P. G. Hoeks, P. J. Brands, F. A. M. Smeets, and R. S. Reneman,“Assessment of the distensibility of superficial arteries,” UltrasoundMed. Biol., vol. 16, no. 2, pp. 121–128, 1990.
[31] H. C. Andrews and B. R. Hunt, Digital Image Restoration. Englewood
Cliffs, NJ: Prentice-Hall, 1977, pp. 148–150.
[32] H. Kanai, H. Satoh, K. Hirose, and N. Chubachi, “A new method formeasuring small local vibrations in the heart using ultrasound,” IEEE
Trans. Biomed. Eng., vol. 40, pp. 1233–1242, 1993.
[33] J. A. Cadzow and O. M. Solomon, Jr., “Linear modeling and thecoherence function,” IEEE Trans. Acoust., Speech, Signal Processing,vol. ASSP-35, pp. 19–28, 1987.
[34] A. A. Luisada and D. M. MacCanon, “The phases of the cardiac cycle,”Amer. Heart J., vol. 83, pp. 705–711, 1972.
[35] Y. Koiwa, R. Hashiguchi, T. Ohyama, S. Isoyama, S. Satoh, H. Suzuki,
and T. Takishima, “Measurement of the instantaneous viscoelastic
properties by impedance frequency curve of the ventricle,” Amer. J.
Physiol., vol. 245, pp. H672–H684, 1986.
[36] R. Hashiguchi, Y. Koiwa, T. Ohyama, T. Takagi, J. Kikuchi, J. P. Butler,
and T. Takishima, “Dependence of instantaneous transfer function onischemic myocardial volume,” Circ. Res., vol. 63, pp. 1003–1011, 1988.

Hiroshi Kanai (A’88–M’91) was born in Matsumoto, Japan, on November 29, 1958. He received
the B.E. degree from Tohoku University, Sendai,
Japan in 1981, and the M.E. and the Dr.Eng. degrees, also from Tohoku University, in 1983 and1986, both in electrical engineering.
From 1986 to 1988, he was with the Education
Center for Information Processing, Tohoku University, as a Research Associate. From 1990 to 1992,
he was a Lecturer in the Department of ElectricalEngineering, Faculty of Engineering, Tohoku University. Since 1992, he has been an Associate Professor in the Department
of Electrical Engineering, Faculty of Engineering, Tohoku University. Hispresent interest is in ultrasonic measurements and digital signal processing
for diagnosis of the heart diseases and arteriosclerosis.
Dr. Kanai is a member of the Acoustical Society of Japan, the Instituteof Electronics Information and Communication Engineering of Japan, theJapan Society of Mechanical Engineers, the Japan Society of Ultrasonics
in Medicine, the Japan Society of Medical Electronics and Biological Engineering, the Institute of Electrical Engineers of Japan, and the Japanese
Circulation Society.

Michie Sato was born in Toyama, Japan, on
November 30, 1972. She received the B.E. degreefrom Tohoku University, Sendai, Japan in 1995 inelectrical engineering. Since 1995, she has been
pursuing the Master’s degree in electrical andcommunication engineering, Tohoku University,
Japan.
Her current interest is in digital signal processing
and ultrasonic measurement for cardiovasculardisease.
Miss Sato is a member of the Acoustical Societyof Japan and the Institute of Electrical Engineers of Japan.

Yoshiro Koiwa was born in Sendai, Japan, in 1944.
He graduated from Tohoku University, Sendai,
Japan, in 1969. He received the M.D. degree fromTohoku University in 1977.
He is presently an Associate Professor of Internal
Medicine at Tohoku University. His main researchinterest is cardiovascular disease, especially cardiac
function and heart failure.
Dr. Koiwa is member of the American Federation
for Clinical Research, the Japanese Circulation
Society, and the Japan Society of Medical
Electronics and Biological Engineering.

Noriyoshi Chubachi (M’83) was born in Kokura,
Japan, on October 5, 1933. He received the B.S.,
M.S., and Ph.D. degrees in electrical engineeringfrom Tohoku University, Sendai, Japan, in 1956,
1962, and 1965, respectively.
In 1965 he joined the Research Institute of Electrical Communication, Tohoku University, where he
was an Associate Professor from 1966 to 1978.
Since 1979, he has been a Professor in the Department of Electrical Engineering, Tohoku University.
From 1982 to 1983, he was a Visiting Professor of
Electrical and Computer Engineering, University of California, Santa Barbara.
He has worked on ultrasonic transducers and delay lines, surface acoustic
devices, acoustoelectronics, piezoelectric materials, acoustic microscopy, andrelated problems.
Dr. Chubachi is a member of the Acoustical Society of America, theInstitute of Electronics and Communication Engineers of Japan, the Society ofJapanese Applied Physics, the Acoustical Society of Japan, the Japan Societyof Ultrasonics in Medicine, the Japan Society for Nondestructive Inspection, the Japan Society of Medical Electronics and Biological Engineering,
the Japan Society of Mechanical Engineers, and the Institute of ElectricalEngineers of Japan. He served as Chairman, Tokyo Chapter of the IEEE
Ultrasonics, Ferroelectrics, and Frequency Control Society from 1987 to 1988.
He is currently serving as Vice President of the Acoustical Society of Japan.