IEICE TRANS. FUNDAMENTALS, VOL.E82–A, NO.4 APRIL 1999

572

PAPER

Special Section on Advanced Signal Processing Techniques for Analysis of Acoustical and Vibrational Signals

PARCORR-Based Time-Dependent AR Spectrum
Estimation of Heart Wall Vibrations
Hiroshi KANAI† , Member and Yoshiro KOIWA†† , Nonmember

SUMMARY
We present a new method for estimation ofspectrum transition of nonstationary signals in cases of low
signal-to-noise ratio (SNR). Instead of the basic functions employed in the previously proposed time-varying autoregressive(AR) modeling, we introduce a spectrum transition constraint
into the cost function described by the partial correlation (PARCORR) coeﬃcients so that the method is applicable to noisy
nonstationary signals of which spectrum transition patterns arecomplex. By applying this method to the analysis of vibration
signals on the interventricular septum (IVS) of the heart, noninvasively measured by the novel method developed in our laboratory using ultrasonics, the spectrum transition pattern is clearlyobtained during one cardiac cycle for normal subjects and a patient with cardiomyopathy.
key words: time-dependent spectrum estimation, AR modeling,
noninvasive diagnosis

1.

Introduction

Much work has been done on parametric spectrum
estimation using the autoregressive (AR) model. A
strong restriction of methods based on this model liesin the necessary assumption that the signals may beconsidered to be stationary during the observation period. Time-varying parametric approaches of modeling have been proposed to overcome this limitationand to take the eﬀects of nonstationary signals into
account explicitly. To estimate the parameters using
a linear algorithm, the unknown time-varying parameters are approximated by linearly weighted combinations of a small number of known functions. The
choice of the basic functions is an important part ofsuch a modeling process. A convenient way is to replace the time-varying coeﬃcients with their secondorder expansion [1] or an arbitrary order expansion [2],[3]. Legendre [4], [5], Fourier [6], prolate spheroidal [7],and B-spline [8] are the basic functions usually chosen.
Since the number of unknown parameters is large, eﬃcient equivalent representations for the modeling, suchas lattice ﬁlters, have also been proposed [2], [7], [9].
However, if the spectrum transition pattern is complex and/or there are large diﬀerences in the transition
Manuscript received August 17, 1998.
Manuscript revised November 11, 1998.
†
The author is with the Department of Electrical Engineering, Graduate School of Engineering, Tohoku University, Sendai-shi, 980–8579 Japan.
††
The author is with the First Department of Internal
Medicine, School of Medicine, Tohoku University, Sendaishi, 980–8578 Japan.

patterns among the individual nonstationary signals, it
is diﬃcult to estimate the transition pattern stably bychoosing a set of basic functions a priori.
We have proposed a method for analyzing the spectrum transition of the multi-frame signals of the fourth
heart sounds detected during the stress test [10]. In thismethod, however, the signals which can be analyzed arelimited to multiple short length signals; the spectrum
transition patterns between these signals are obtained.
In this paper, by modifying this method, we propose
a new approach of modeling to estimate the spectrum
transition of a nonstationary signal by using a linear
algorithm without any basic functions.
In this paper, moreover, we describe the spectrum
transition constraint not by the linear predictive coefﬁcients of the AR model but by the partial correlation(PARCORR) coeﬃcients. In the method developed in[10], the cost function of the multi-frame signals is deﬁned by the sum of the residual powers in the AR modeling and the spectrum transition between the multiframe signals. The latter component of the spectrum
transition is deﬁned by the sum of the diﬀerences ofthe linear predictive coeﬃcients {ai } between the successive frames. In general, however, the values of thecoeﬃcients {ai } are large at low order and small at highorder. Thus, it is signiﬁcant to determine the weightof each order of the coeﬃcients {ai } in the summation
process.
For this problem, in this paper, the PARCORR
coeﬃcients {ki } instead of the linear predictive coeﬃcients {ai } are introduced into the deﬁnition of the costfunction. Since the range of the PARCORR coeﬃcients{ki } is from −1 to +1 even in low and high orders, it is
not necessary to consider the weights in the deﬁnitionof the spectrum transition of the cost function.
To noninvasively diagnose the acoustic characteristics of the heart muscle, it is necessary to measure thesmall vibration signals on the heart wall from the chest
surface with enough accuracy to permit analysis in thefrequency range up to at least 100 Hz and analysis of
the resultant nonstationary signal continuously during
one cardiac cycle.
For the problem of the noninvasive measurement
described above, we have already developed a novel
method for accurate ultrasonic-based measurement ofsmall velocity signals on and in the heart wall in thefrequency range up to several hundred Hertz [11]–[13].

KANAI and KOIWA: SPECTRUM ESTIMATION OF HEART WALL VIBRATIONS

573

If non-stationary spectrum analysis is applied to
the resultant velocity signals, more detailed information on the heart function can be obtained. To solve
the problem of the spectrum analysis, we apply thenovel time-varying modeling developed in this paper
to the nonstationary small vibration signals on the interventricular septum (IVS) to diagnose the acoustic
characteristics of the heart muscle. These characteristics and the transition patterns show the potential of
the proposed method for the noninvasive diagnosis ofheart diseases.
2.

Principle of Spectrum Estimation Using
PARCORR Coeﬃcients

Let us divide an original nonstationary signal x(n) into
successive F short signals {xj (n)}, n = 0, 1, . . . , N − 1,
j = 0, 1, . . . , F − 1, each called a frame, where F isthe number of frames. Let us assume that each frame
signal xj (n) is an AR signal of order M , represented
by the forward and backward recursions. The order Mof the AR model coincides with double the number ofthe poles in the frequency range up to the Nyquist frequency. The number of poles is physically determinedby the eigen-vibration systems. During the whole analyzed period of the original nonstationary signal x(n),let us assume that the eigen-vibration systems do notvary and that only their eigen-frequencies vary. In thispaper, therefore, we assume that the all of the frame
signals are of the same order M .
The forward predictive error efM,j (n) and the backward predictive error ebM,j (n) of M -order prediction arerespectively given by
efM,j (n) =MaMi,j · xj (n − i),

i=0
M


bMi,j · xj (n − i − 1),(2)

where {aMi,j } and {bMi,j } are ith order forward andbackward predictive coeﬃcients, respectively, of jthframe data, aM0,j = 1 and bMM,j = 1. When the predictive order is equal to M , the average power αM,j ofthe predictive error for the data in the period [M, N −1]
of jth frame data is given by(i = 1, 2, . . . , M )Equation (5) is rewritten as follows:MCi ,j · aM ,j = −Ci0,j . (i = 1, 2, . . . , M )

(6)

=1

For the backward prediction, on the other hand, it
is well-known [14], [15] that the backward predictive error signal ebM−1,j (n− M ) in the (M − 1)th order prediction should be orthogonal to the (M − 1)th order back
ward predictor, xbM−1,j (n−M ) = M−2i=0 bM−1,i,j xj (n−1 − i), of the jth frame signal, that is,

N
−1

1
ebM−1,j (n − M )
N −M
n=M
M−2


×
bM−1,i,j xj (n − 1 − i)
= 0.
(7)
i=0

From this equation, for the cases of p > q, the error
series ebp−1,j (n − p) is orthogonal to each of the components xj (n − p + 1), xj (n − p + 2), . . . , xj (n − 1) of the
backward predictor x
bp−1,j (n − p). Thus,
N −1


1
ebp−1,j (n − p) · ebq−1,j (n − q) = 0.
N −M

 f
1
|eM,j (n)|2
N −M
n=M

M 
M


aMi,j aM ,j Ci ,j ,(3)

i=0 =0where Ci ,j is the covariance matrix of data xj (n) deﬁned by
1
N −M

N
−1

n=M

xj (n − i) · xj (n − ).

(p > q) (8)For the case of p = q, however, by deﬁning a nonnegative constant δp,j , the following relation is obtained:δp,j =

N
−1

 b

1
ep−1,j (n − p)2
N −M
n=M

=

p−1
p−1


bp−1,i,j bp−1, ,j Ci+1, +1,j .

(9)i=0 =0Moreover, it is well-known [14], [15] that the forward predictive coeﬃcient aMi,j is recursively derived
from the results estimated up to the (M − 1)th order
as follows:

N −1

Ci ,j =

(5)

=0

(1)

i=0

=

M


1 ∂αM,j=
aM ,j · Ci ,j = 0
2 ∂aMi,j

n=M

ebM,j (n − M − 1) =αM,j =

By minimizing αM,j with respect to the forward predictive coeﬃcients {aMi,j }, the resultant normal equation
is given by(4)

aMi,j = aM−1,i,j + kM,j bM−1,i−1,j ,(i = 1, 2, . . . , M − 1)

(10)where the term kM,j is the PARCORR coeﬃcient forjth frame data. From this recursion,aMi,j = (aM−2,i,j + kM−1,j bM−2,i−1,,j )IEICE TRANS. FUNDAMENTALS, VOL.E82–A, NO.4 APRIL 1999

574

+ kM,j bM−1,i−1,j= ···= (a1i,j + k2,j b1,i−1,,j ) + k3,j b2,i−1,,j + · · ·+ kM,j bM−1,i−1,,j= k1,j +M

k ,j b −1,i−1,j ,

(11)

=2

where the relation a1i,j = k1,j is employed.
Let us respectively deﬁne M -dimension vectors aj ,kj , cj , an (N − M )-dimension vector xj , an M × M
upper triangle matrix Bj , an M × M diagonal matrix∆j , and an M × M covariance matrix Cj , which is
positive deﬁnite, by
aj = [aM1,j , aM2,j , · · · , aMM,j ]T ,T

kj = [kM1,j , kM2,j , · · · , kMM,j ] ,

(12)
(13)

T

cj = [C10,j , C20,j , · · · , CM0,j ] ,(14)1
xj = √[xj (M ), · · · , xj (N − 1)]T ,
(15)
N −M


1 b10,j b20,j . . .
bM−1,0,j
 0
1
b
.
.
.
bM−1,1,j 
21,j




..
.
.

, (16).
0
1.
Bj =  0

 ..

..
..
..
 .
.
.
. bM−1,M −2,j 
0 ...
...
0
1


0
δ1j 0 . . .
 0 δ2j . . .
0 ∆j =  .
(17)
.. ,
.
.
 ..
.
0
. 
0 . . . 0 δM,j
Cj = [Ci ,j ].

(18)

Using these vectors and matrices, Eqs. (6), (11), and(8), Eq. (9) are respectively simpliﬁed as
Cj · aj = −cj ,aj = Bj · kj ,

(19)
(20)

BTj · Cj · Bj = ∆j .

(21)

Thus, the total power αM,j of the M th-order forward
predictive error in Eq. (3) is given by using the PARCORR coeﬃcients as follows:
αM,j = aTj Cj aj + 2cTj aj + xTj xj= kTj BTj Cj Bj kj + 2cTj Bj kj + xTj xj= kTj ∆j kj + 2(BTj cj )T kj + xTj xj .
3.

(22)

Principle of Estimation of Spectrum Transition

When each frame of the multi-frame nonstationary signal {xj } (j = 0, 1, . . . , F − 1) is described by the AR
model, based on Eq. (22), let us deﬁne the logarithmic
likelihood function , which shows the probability of{xj } for unknown PARCORR coeﬃcients {kj }, byF−1kTj ∆j kj + 2(BTj cj )T kj + xTj xj=−|xj |2j=0

+ λj |kj+1 − kj |2 ,(23)

where {λj } are F Lagrange multipliers. In Eq. (23),
we assume that the frame data are cyclic and then 0th
frame data exactly coincide with the F th frame data,that is, x0 (n) = xF (n) and k0 = kF . The second
term in the right-hand side of Eq. (23) shows the constraint for the spectrum transition between the successive frames. The denominator |xj |2 of the ﬁrst term in
the right-hand side normalizes the signal power of jthframe data.
In this logarithmic likelihood function , if {λj } are
set to be large, the second term of Eq. (23) is dominant
and then the maximization of  corresponds to the solu2
tion so as to be |kj+1 − kj | = 0, for j = 0, 1, . . . , F −1.
However, this procedure does not at all take into consideration the spectrum estimation which is achievedby the minimization of αM,j in Eq. (22).
On the other hand, if the values of {λj } are set tobe very small, the second term of Eq. (23) is negligible
and then only its ﬁrst term is dominant. Thus, it isimportant to select appropriate values of {λj }. In thispaper, the values are set uniformly for the subjects in invivo experiments so that the resultant spectrum transition patterns become reasonable from the cardiologicalside and the diﬀerences between the healthy subjects
and the patient are clear.
Let us determine the PARCORR coeﬃcients {kj }which maximize the logarithmic likelihood function 
as follows. Let us assume that F Lagrange multipliers{λj } are the same value of λ. By taking the partial
derivative of  with respect to {kj } and setting theresults to be zero, the following simultaneous equations
are obtained.
1 ∂
−
2 ∂kj

1 
∆j kj + BTj cj + λ (2kj − kj−1 − kj+1 )
=
2
|xj |
= 0.
(24)
By solving Eq. (24), the PARCORR coeﬃcients {kj }of all frame data are estimated under the constraint forthe spectrum transition between the frame data.
From Eq. (24),
BTj cj
∆j
+2λI
k
−λk
−λk
=
−
, (25)jj−1
j+1|xj |2|xj |2where I is the M × M identity matrix. Let us denote(∆j /|xj |2 + 2λI) by an M × M diagonal matrix Dj
and let us deﬁne an F M × F M matrix G and F M dimension vectors g and k as follows:KANAI and KOIWA: SPECTRUM ESTIMATION OF HEART WALL VIBRATIONS

575



D0 −λI
0···0 −λI D1 −λI 0
·
··


.
.
..
..
 0
−λI D2

G= .
.
.
.
..
..
..
 ..
0


..
..
..
..
 0
.
.
.
.
−λI
0
...
0 −λI

−λI
0
..
.






,

0


−λI 
DF −1


T  T T T
T
B
c
B
c
F
−1
0
F −1
0
 ,
g=
,···,
|xj |2|xF −1 |2T
k = kT0 , kT1 , . . . , kTF −1 .



(26)
(27)
(28)

Using these matrix and vectors, optimum estimation of
the spectrum transition of a nonstationary signal under the constraint for the spectrum transition between
frame data is simultaneously achieved by solving thefollowing linear simultaneous equations:
Gk = −g,(29)

where the term −λI in the top right and bottom leftof the matrix G is introduced by assuming x0 (n) =xF (n). Therefore, the PARCORR coeﬃcients of themulti-frame data are simultaneously estimated by = −G−1 g,k(30)

where G−1 is the inverse matrix of G. Using the estiTT , . . . , kT
= kT , k, estimates of the vecmates k
0

1

F −1

tors aj of linear predictive coeﬃcients in the jth frame
data are obtained by
j . (j = 0, 1, . . . , F − 1)aj = Bj k
4.

(31)

In vivo Experiments for Estimation of theSpectrum Transition of the Vibrations on
the Heart Wall

We applied this method to the analysis of the velocity
signals measured in our laboratory using a noninvasive
measurement method [12] on the interventricular septum (IVS) in the heart wall of two healthy male subjects and a patient with acute lymphoblastic leukemia
and serious doxorubicin-cardiomyopathy, who had beentreated with an anti-cancer drug (antracenadiones).
Figure 1 shows a standard brightness (B)-mode
short-axis image of the cross-sectional area around thedetected points preset in the IVS of a presumedly
healthy 26-year-old male volunteer. The ultrasonicbeam passing through the two points of R and L isalmost perpendicular to the IVS during the measurements.

Fig. 1 A standard B-mode short-axis image showing the crosssectional area around the detected points preset in the interventricular septum (IVS) in an in vivo experiment for the detectionof their instantaneous positions and the velocity signals of a presumedly healthy 26-year-old male volunteer. Points (R) and (L)are on the RV surface and the LV surface of the IVS, respectively.

4.1 Waveforms for Spectrum Analysis
Figures 2 (1-a) and (1-b) show the electrocardiogram(ECG) and the heart sound (PCG), respectively, of the
26-year-old healthy male subject in Fig. 1. Figure 2 (1c) shows the velocity signals {vi (t)} (i = 0, 1, . . . , 13)
measured for the 14 points preset from the right ventricular (RV) surface (the point R in Fig. 1) to the left
ventricular (LV) surface (the point L in Fig. 1) of the
IVS. The 14 velocity signals are overlaid.
From the signals in Fig. 2 (1-c), the time
diﬀerential waveforms {vi (t)} of the velocity signals{vi (t)}, which are the acceleration signals, are obtained
to emphasize higher frequency components as shown in
Fig. 2 (1-d). The original sampling frequency is 4.5 kHz.
The signals {vi (t)} are re-sampled at a frequency of
200 Hz so that they are analyzed in the frequency range
up to 100 Hz. Figure 2 (2) shows the waveforms for the
32-year-old male patient.
One heartbeat signal vi (t) on the LV side of the
IVS in Figs. 2 (1-d) and (2-d) is divided into multi
ple frame signals {vi (n; j)}, (j = 0, 1, . . . , F − 1; n =0, 1, . . . , N −1), each of which is N = 20 points (100 ms)
in length, which roughly corresponds to the duration of
the ﬁrst heart sound. Adjacent frame signals overlap
each other by a three-quarter length (75 ms). For thesignals in Figs. 2 (1-d) and (2-d), the numbers (F ) of
the preset frames are 47 and 33, respectively, which are
determined by the length of one cardiac cycle.
4.2 Estimates of the Spectrum Transition
By applying the proposed method (M = 6) to the accel-

IEICE TRANS. FUNDAMENTALS, VOL.E82–A, NO.4 APRIL 1999

576

tively. These results correspond to those obtained byindependently applying the standard AR estimation
method to each frame signal. By comparing the estimates in Figs. 3 (a) and (c) with those in Figs. 4 (a)and (c), the introduction of the constraint in Eq. (23)
is found to be eﬀective.
4.3 Discussion

Fig. 2 (a) ECG, (b) PCG, (c) velocity signals {vi (t)}, (d) ac
celeration signals {vi (t)} of (c). (1) for the healthy male subject
in Fig. 1 and (2) for a 32-year-old male patient.

eration signals v13 (t) on the LV surface in Fig. 2 (1-d),
the resultant spectrum transition patterns are shown in
Fig. 3 (a). The center frequency of each estimated pole
is indicated by a circle. In these experiments, the valueof λ of Eq. (23) is set as 0.316.
Figure 3 (b) shows the results obtained by applying
the same procedure to the acceleration signal of a 22year-old healthy male subject.
On the other hand, Figs. 3 (c) and (d) show the
spectrum transition patterns of the acceleration signal
on the LV surface of the IVS for a 32-year-old malepatient. The data analyzed in Figs. 3 (c) and (d) (also in
Fig. 2 (2)) were measured ﬁve months and two months
before his death, respectively.
For the same acceleration signals in Figs. 3 (a) and(c), the instantaneous spectrum patterns are estimatedwithout using the transition constraint, that is, λ = 0.
The results are shown in Figs. 4 (a) and (c), respec-

1. For the data in Figs. 3 (c) and (d) of the patient,
it is found that the duration of diastole becomesshorter, while that of systole is still almost thesame as those of the normal subjects in Figs. 3 (a)and (b).
2. Since clear components of the heart sound havenot been recognized for healthy subjects in the diastole from the end of the second heart sound (II)to the beginning of the ﬁrst heart sound (I), analysis of the heart sound in diastole has not beenreported in the literature. As shown by the ovals
in Figs. 3 (a)–(d), however, there are the remarkable spectrum peaks between 10 Hz and 20 Hz forboth the normal subjects and the patient. Moreover, the dominant frequencies gradually become
lower from the beginning to the end of diastole,
which corresponds to the dilation of the LV.
3. By comparing the results in Figs. 3 (a) and (b) withthose in Figs. 3 (c) and (d), the power of the normal subjects is seen to be large in the frequency
range up to 100 Hz around the radiation timing of
the second heart sound (II), which corresponds tothe period from the end-systole to the beginningof diastole.
For the patient, however, the power decreased,
especially in the higher frequency band. Thesephenomena correspond to the myocardial defectwhich occurred due to the dose of the anti-cancer
drug [13].
4. In systole, since there is a large diﬀerence
among individual subjects, remarkable characteristics have not yet been recognized. For the radiation timing of the ﬁrst heart sound (I), however,the power of the patient in Figs. 3 (c) and (d) decreases, especially at higher frequencies.
5.

In vivo Experiments for Estimation of theSpatial-Transition of the Spectra of the Vibrations on the Heart Wall

Figures 5 (1) and (2) show the acceleration signals at
the 14 points preset at even intervals of 0.75 mm fromthe RV side to the LV side of the IVS for the normal
subject in Fig. 2 (1) and the patient in Fig. 2 (2), respectively. The one cardiac signal is divided into several
periods, each of which is shown in these ﬁgures.
Figures 6 (1) and (2) show the spectrum transition

KANAI and KOIWA: SPECTRUM ESTIMATION OF HEART WALL VIBRATIONS

577

Fig. 3 Spectrum transition pattern of the acceleration signal on the heart wall for one
cardiac cycle (circle: the center frequency of each estimated pole). (a) the 26-year-old
healthy male subject, (b) a 22-year-old healthy male subject, (c) a 32-year-old male patient(ﬁve months before death), (d) the same 32-year-old male patient (two months before
death). “I” and “II” show the radiated timing of the ﬁrst and the second heart sounds.

Fig. 4 Spectrum transition pattern of the acceleration signal
estimated without using the transition constraint, that is, λ = 0.
(a) the 26-year-old healthy male subject in Fig. 3 (a), (c) the 32year-old male patient (ﬁve months before death) in Fig. 3 (c).

Fig. 5 The acceleration signals from the RV to the LV and theiranalyzed periods from A to E or from A to D. Upper: for thenormal subject in Fig. 2 (1). Lower: for the patient in Fig. 2 (2).

IEICE TRANS. FUNDAMENTALS, VOL.E82–A, NO.4 APRIL 1999

578

Fig. 6 The spatial transition pattern of the estimated spectrafrom the RV to the LV for each period from A to E or from A
to D indicated in Fig. 5 (circles: the center frequency of eachestimated pole). Upper 5 ﬁgures: for the normal subject. Lower
4 ﬁgures: for the patient.

patterns obtained by applying the proposed method(M = 6) to the sectional signals in each of the periods from A to E of Fig. 5 (1) and from A to D of 5 (2),respectively.
For the 3 subjects, the average value ∆Pi of thediﬀerence of the power spectrum Pi (f ) of the ith layerfrom the power spectrum PR (f ) of the RV side in thefrequency range from d.c. to 100 Hz is evaluated for theperiods of A and C as follows:
1 (Pi (f ) − PR (f )).
101
100

∆Pi =

(32)

f =0

The results are shown in Fig. 7.
From Figs. 6 and 7, for the normal subject, there
are large changes in power among the signals from theRV side to the LV side, especially in period A (at the beginning of the systole), period C (from the end-systole
to the beginning of the diastole), and period D (diastole). That is, the power on the LV side is larger thanthe power on the RV side. This increase in power showsthat there is a large change in thickness in the IVS during each period for the normal subject.

Fig. 7 The average value ∆Pi of the diﬀerence of the power
spectrum Pi (f ) of the ith layer from the power spectrum PR (f )of the RV side in the frequency range from d.c. to 100 Hz for four
subjects in Fig. 3. upper: for the period A around the R-wave.
lower: for the period C around the second heart sound (II).

For the patient, however, such remarkable phenomena were not observed in any period, which corresponds to the decrease in thickness of the IVS in theserious patient.
If the myocardium is a passive component, thesevalues of the change in thickness directly correspondto the strain in the heart wall and its elasticity may
be evaluated using the blood pressure. However, since
the myocardium itself is not a passive component, butrather an active component, it will be diﬃcult to estimate the elasticity of the myocardium from the obtained change in thickness in this paper.
6.

Conclusions

We have presented a new method for estimation of thespectrum transition of a nonstationary signal in low
SNR cases using a linear algorithm. Applying the proposed method to heart wall vibrations, we found clear
spectrum transition patterns.
The small velocity signals accurately measured by

KANAI and KOIWA: SPECTRUM ESTIMATION OF HEART WALL VIBRATIONS

579

our method contain suﬃcient information to diagnose
the acoustic characteristics of the heart muscle during
one cardiac cycle. Thus, the measurement of the heart
wall vibrations and their analysis as proposed in thispaper are expected to lead to the development of a newscientiﬁc ﬁeld of noninvasive diagnosis of heart dysfunction.
In this paper, the order M of the AR model was
ﬁxed to be 6 based on our previous analysis of the
fourth heart sound in [10], and λ = 0.316 was employed
as the value of the Lagrange multipliers based on invivo experiments using two healthy subjects and one
patient. It is, however, necessary to develop a methodfor deciding the order M of the AR model and valuesof the Lagrange multipliers. These problems are underinvestigation.

method for measuring small local vibrations in the heart
using ultrasound,” IEEE Trans. BE, vol.40, no.12, pp.1233–
1242, 1993.
[12] H. Kanai, M. Sato, Y. Koiwa, and N. Chubachi, “Transcutaneous measurement and spectrum analysis of heart wall
vibrations,” IEEE Trans. UFFC, vol.43, no.5, pp.791–810,
1996.
[13] H. Kanai, H. Hasegawa, N. Chubachi, Y. Koiwa, and M.
Tanaka, “Noninvasive evaluation of local myocardial thickening and its color-coded imaging,” IEEE Trans. UFFC,
vol.44, no.4, pp.752–768, 1997.
[14] J.D. Markel and A.H. Gray, Jr., “Linear prediction of
speech,” Springer-Verlag, 1978.
[15] S.M. Kay, “Modern spectral estimation,” Prentice Hall,
Chapter 6, 1987.

Acknowledgments
The authors would like to thank Prof. Noriyoshi
Chubachi of Tohoku Gakuin University, Prof. Motonao
Tanaka of Tohoku Kousei-Nenkin Hospital, Dr. Yoshiko
Saito of Tohoku University School of Medicine, Dr.
Yoshiro Uzuka of the Sendai Center for Hematologic
Disorders, and Miss Michie Sato for their helpful comments and their support in the in vivo experiments.
References[1] T.S. Rao, “The ﬁtting of non-stationary time-series models
with time-dependent parameters,” J. of the Royal Statist.
Soc. Series B, vol.32, no.2, pp.312–322, 1970.
[2] L.A. Liporace, “Linear estimation of non-stationary signals,” J. Acoust. Soc. Amer., vol.58, no.6, pp.1288–1295,
1975.
[3] K.C. Sharman and B. Friedlander, “Time-varying autoregressive modeling of a class of nonstationary signals,” Proc.
IEEE ICASSP, pp.22.2.1–22.2.4, 1984.
[4] F. Kozin, “Estimation and modeling of nonstationary timeseries,” Proc. Symp. Appl. Comput. Meth. Eng., vol.1, Los
Angeles, CA, pp.603–612, 1977.
[5] G. Alehgrin, M. Barlaud, and J. Menez, “Unbiased parameter estimation of nonstationary signals and noise,” IEEE
Trans. ASSP, vol.ASSP-34, no.5, pp.1319–1322, 1986.
[6] M. Hall, A.V. Oppenheim, and A. Willsky, “Time-varying
parametric modeling of speech,” Proc. IEEE Decision andControl Conf., New Orleans, LA, pp.1085–1091, 1977.
[7] Y. Grenier, “Time-dependent ARMA modeling of nonstationary signals,” IEEE Trans. ASSP, vol.ASSP-31, no.4,
pp.899–911, 1983.
[8] M.J. Flaherty, “Application of polynomial splines to atime-varying autoregressive model of speech,” Proc. IEEEICASSP, pp.2220–2223, 1988.
[9] E. Karlsson and M.H. Hayes, “ARMA modeling of timevarying systems with lattice ﬁlters,” Proc. IEEE ICASSP,
pp.2335–2338, 1986.
[10] H. Kanai, N. Chubachi, K. Kido, Y. Koiwa, T. Takagi,
J. Kikuchi, and T. Takishima, “A new approach to timedependent AR modeling of signals and its application to
analysis of the fourth heart sound,” IEEE Trans. SP, vol.40,
no.5, pp.1198–1205, 1992.
[11] H. Kanai, K. Hirose, H. Sato, and N. Chubachi, “A new

Hiroshi Kanai
was born in Matsumoto, Japan, on November 29, 1958. Hereceived the B.E. degree from Tohoku
University, Sendai, Japan in 1981, and the
M.E. and the Dr.Eng. degrees, also fromTohoku University, in 1983 and 1986,
both in Electrical Engineering. From1986 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 Electrical Engineering, Faculty ofEngineering, Tohoku University. Since 1992, he has been an Associate Professor in the Department of Electrical Engineering,
Faculty of Engineering, Tohoku University. His present interest
is in ultrasonic measurements and digital signal processing for
diagnosis of the heart diseases and atherosclerosis. Dr. Kanai ismember of the Acoustical Society of Japan, the Japan Society of
Ultrasonics in Medicine, the Japan Society of Medical Electronics
and Biological Engineering, the Institute of Electrical Engineersof Japan, the Japanese Circulation Society, the Japanese College
of Cardiology, and the IEEE.

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