Fast characterization of two ultrasound longitudinal waves
in cancellous bone using an adaptive beamforming technique
Hirofumi Takia)
Graduate School of Engineering, Tohoku University, Sendai 980-8579, Japan

Yoshiki Nagatani
Department of Electronics, Kobe City College of Technology, Kobe 651-2194, Japan

Mami Matsukawa
Faculty of Science and Engineering, Doshisha University, Kyotanabe 610-0321, Japan

Katsunori Mizuno
Institute of Industrial Science, The University of Tokyo, Tokyo 113-8654, Japan

Toru Sato
Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan(Received 6 September 2014; revised 18 January 2015; accepted 4 March 2015)
The received signal in through-transmission ultrasound measurements of cancellous bone consists of
two longitudinal waves, called the fast and slow waves. Analysis of these fast and slow waves may
reveal characteristics of the cancellous bone that would be good indicators of osteoporosis. Because
the two waves often overlap, decomposition of the received signal is an important problem in thecharacterization of bone quality. This study proposes a fast and accurate decomposition method
based on the frequency domain interferometry imaging method with a modified wave transfer function that uses a phase rotation parameter. The proposed method accurately characterized the fast andslow waves in the experimental study, and the residual intensity, which was normalized with respect
to the received signal intensity, was less than 20 dB over the bone specimen thickness range from
6 to 15 mm. In the simulation study, the residual intensity was less than 20 dB over the specimen
thickness range from 3 to 8 mm. Decomposition of a single received signal takes only 5 s using a laptop personal computer with a single central processing unit. The proposed method has great potential
to provide accurate and rapid measurements of indicators of osteoporosis in cancellous bone.
C 2015 Acoustical Society of America. [http://dx.doi.org/10.1121/1.4916276]
V[GH]

Pages: 1683–1692

I. INTRODUCTION

Osteoporosis is a disease that is characterized by adecrease in bone mass, leading to an increased possibility of
the risk of bone fracture. Dual x-ray absorptiometry is thegold standard for osteoporosis diagnosis; however, diagnosis
using this method is based on only one parameter related tobone porosity (Genant and Jiang, 2006; Blake andFogelman, 2007). Additional testing using another modality
may improve the accuracy of the diagnosis of osteoporosis.
Quantitative ultrasound (QUS) is one technique thatoffers the potential to estimate the possibility of osteoporosis
without using ionizing radiation (Langton et al., 1984;
Laugier, 2008; Barkmann et al., 2008a,b). QUS measures
the speed of sound and the broadband ultrasonic attenuation
characteristics that provide important information for characterization of bone quality and strength (Bouxsein et al.,
1995; Njeh et al., 1997; Han et al., 1997; Lochm€uller et al.,
1999; Ha€ıat et al., 2009). It has been established that thereceived signal in through-transmission ultrasound measurements of cancellous bone consists of the fast and slow waves(Hosokawa and Otani, 1997, 1998; Hughes et al., 1999;a)Author to whom correspondence should be addressed. Electronic mail:taki@ecei.tohoku.ac.jp

J. Acoust. Soc. Am. 137 (4), April 2015

Kaczmarek et al., 2002; Cardoso et al., 2003; Fellah et al.,2004). Several studies have reported that analysis of thesefast and slow waves may reveal the characteristics of cancellous bone (Cardoso et al., 2003; Ha€ıat et al., 2007; Hughes
et al., 2007; Hosokawa, 2009; Mizuno et al., 2008; Mizuno
et al., 2009). Because the characteristics of cancellous boneprovide good indications of the bone changes caused byosteoporosis (Njeh et al., 1999), fast and accurate characterization of the fast and slow waves may be important to therealization of an osteoporosis screening method.
Analysis of the fast and slow waves requires decomposition of the received signal from cancellous bone into the fastand slow waves. Because the two waves often overlap,
decomposition of the received signal is an important problem
in the characterization of bone quality using QUS. One proposed strategy is the application of a Bayesian method using a
Markov chain Monte Carlo method with simulated annealing(Marutyan et al., 2007; Anderson et al., 2010; Nelson et al.,2011; Hoffman et al., 2012). The Bayesian decomposition
method succeeded in accurately characterizing the fast andslow waves of a received signal from cancellous bone; however, this method involves a considerable computational load.
Other groups have reported fast decomposition methods
for characterization of the fast and slow waves. The spacealternating generalized expectation-maximization method0001-4966/2015/137(4)/1683/10/$30.00

C 2015 Acoustical Society of America
V

1683

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21

has been applied to signals that simulated the received signal
from a proximal femur (Dencks and Schmitz, 2013). Wearreported that both the modified least squares Prony’s methodwith curve fitting (Wear, 2013) and band-limited deconvolution (Wear, 2014) are suitable for rapid characterization of
the fast and slow waves. However, these fast decomposition
methods have only modest accuracy when compared withthat of the Bayesian decomposition method.
In the present study, we propose a fast decomposition
method based on frequency domain interferometry (FDI),
with decomposition accuracy that is similar to that of a
Bayesian method. We explain the methodology of the proposed decomposition method, present experimental results,
and offer discussion and conclusions related to our findings.
II. METHODS

The proposed decomposition method requires no priorinformation about the characteristics of a bone specimen,e.g., the specimen thickness. The method uses a specific parameter set for the bone characteristics, and estimates the
specimen thickness as part of the decomposition process.
After decomposition, we can estimate the characteristics ofthe bone specimen using its true thickness and the estimatedthickness. In this section, we explain the modified transfer
functions used for the fast and slow waves that satisfy the
Kramers–Kronig relations. We then describe the proposed
decomposition method based on the FDI imaging method.
A. Transfer functions for fast and slow waves

The wave propagation in cancellous bone has been modeled in previous studies (Marutyan et al., 2006; Anderson
et al., 2008). The conventional propagation models are
expressed using the following formula:SR ðf Þ ¼ SI ðf Þ½H1 ðf Þ þ H2 ðf Þ;(1)

where SR(f) is a frequency component of the received signal
passing through a bone specimen in water at frequency f,
SI(f) is the same component passing through a water-only
path, and H1(f) and H2(f) are the transfer functions for the
fast and slow waves, respectively. Several researchers haveused a transfer model that assumes that the fast and slow
waves are propagating simultaneously through a medium
exhibiting linear-with-frequency attenuation (Nelson et al.,2011; Wear, 2013, 2014). This model introduces the assumption that the wavefront of an ultrasound wave passing
through a bone specimen in water is the same as that passing
through a water-only path. However, the heterogeneity ofcancellous bone causes complicated changes in the wavefronts as they pass. Because the ultrasound waves are typically received by a flat transducer, an uneven wavefront has
an effect on the transfer function. In this study, we approximate this effect using the following formula:ð0sR ðtÞ ¼ ½fT1 ðsÞsR1 ðt  sÞ þ fT2 ðsÞsR2 ðt  sÞds;(2)
where sR0 (t) is the received signal of the proposed model inthe time domain, sR1(t) and sR2(t) are the fast and slow wave
1684

J. Acoust. Soc. Am., Vol. 137, No. 4, April 2015

signals of the conventional model in the time domain,
respectively, and fT1(t) and fT2(t) are the functions that
express the arrival time distributions of the wavefront at thereceive element for the fast and slow waves, respectively.
The frequency component of the received signal in the proposed method is expressed as
SR 0 ðf Þ ¼ F½sR 0 ðtÞ ¼ SI ðf Þ½FT1 ðf ÞH1 ðf Þ þ FT2 ðf ÞH2 ðf Þ;(3)
where F denotes a Fourier transform, and FT1(f) and FT2(f)are the frequency components of the arrival time distributions of the fast and slow waves, respectively. We call thiseffect, which is caused by the arrival time distribution of the
wavefront, the path-length-variation (PLV) effect.
Figure 1 shows an example of the arrival time distribution of the wavefront at the receive element, and of the PLV
effect on the intensity and phase in the frequency domain. In
this case, the distribution follows a chi-squared distribution
with three degrees of freedom. Typically, the phase of FTi(f)is not proportional to f and the intensity of FTi(f) is not constant, where i ¼ 1 and 2 for the fast and slow waves, respectively. Therefore, we approximate FTi(f) in the mainfrequency band using the following formula:
FTi ðf Þ ﬃ ATi exp ½ci f þ jðdi f þ hi Þ;(4)where ATi, ci, ri, and hi are real constants. cif and jrif denotethe attenuation and the time shift caused by the arrival time
distribution at the receive element, respectively. We call hi
the phase rotation parameter, and it is independent of f. The
modified transfer function with the additional phase rotation
parameter is expressed by the following formula:
Hi 0 ð f Þ ¼ FTi ð f ÞHi ð f Þ
2pfd 2pfd
¼ Ai 0 exp bi 0 fd þ jþ jdi f þ jhi ;ci ð f Þ cW
(5)
Ai 0 ¼ Ai ATi ;(6)cbi 0 ¼ bi þ i ;d(7)where Ai is the signal amplitude parameter that is independent of f, bi is the slope of attenuation, i.e., broadband ultrasonic attenuation parameter, in the conventional model, d isthe bone specimen thickness, ci(f) is the phase velocity, andcW is the sound velocity in water. Ai0 and bi0 are real constants, and jdif denotes the time shift. Therefore, the major
modification made in the proposed model is the introductionof the phase rotation parameter hi. The introduction of thisparameter compensates for the waveform change caused by
the uneven wavefront.
In this study, we assume that cW has a constant value.
When the medium has linear-with-frequency attenuation, thephase velocities of the fast and slow waves can be expressedby the following formula (Waters et al., 2005):
 
1
1
bi
f
¼  2 ln;(8)pc i ð f Þ c i ðf 0 Þ
f0Taki et al.: Fast characterization of bone wave

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21

FIG. 1. (a) Example of the arrival time
distribution of the wavefront at thereceive element, and (b) its intensityand (c) phase in the frequency domain.
In this case, the distribution follows achi-squared distribution with threedegrees of freedom. We employ linearwith-frequency models to approximate
the effect of the arrival-time-distribution on the intensity and phase rotationin the major frequency band. ATi and hi
are the intensity and phase rotation parameters given in Eq. (4), respectively.

where f0 is a reference frequency.
We applied the proposed method to simulated and experimental data measured in a previous work (Nagatani et al.,2008). To prepare SI(f), we used the received signal passing
through water in an acoustic tube, where the center frequencyof the pulse was 1 MHz and the pulse length was equal to one
wavelength at the center frequency. Planar polyvinylidene fluoride transducers were used for both transmitter and receiver.
In this experiment, the bone specimen was immersed in waterin the acoustic tube. In the simulation study, the bone specimen size was 17  17  d mm, where the bone specimen
thickness d varied from 2 to 9 mm. In the experimental study,
the bone specimen size was 20  20  d mm, where the bone
specimen thickness d varied from 6 to 15 mm.
The proposed decomposition method used propagation
parameters that were reported in a previous study (Nelson
et al., 2011) as the provisional parameter set. This study also
used the experimental data reported by Nagatani et al.
(2008), and the parameters were estimated by application ofa Bayesian decomposition method to the received signal fora bone specimen thickness of 9 mm. This setting modifies
the transfer function given by Eq. (5) to become
2pfdTi 2pfdTi
Hi 0 ð f Þ ﬃ Ai 0 exp bPi fdTi þ j
cPi ð f ÞcW
(9)þjdi 0 f þ jhi 0 ;
 
1
1
bPi
f
;

¼  2 lnpcPi ð f Þ ci ðf0 Þ
f0

(10)

where dTi is the temporal specimen thickness, and bPif andcPi(f) denote the provisional attenuation coefficient and thephase velocity that were reported in Nelson et al. (2011),
J. Acoust. Soc. Am., Vol. 137, No. 4, April 2015

respectively. In this study, bP1 and bP2 were 49.2 and 7.1 dB/
cm/MHz, respectively, f0 was 1 MHz, and cP1(f0) and cP2(f0)were 1933 and 1475 m/s, respectively. The proposed decomposition method does not require A1 and A2, because it estimated the parameters in a decomposition process. The
temporal specimen thickness dTi is never the same as theestimated thickness, as emphasized in this study, because theoptimum value of dTi satisfies the equation bPidTi ¼ bi0 d.
That is, the use of dTi can compensate for the difference
between the optimum attenuation bi0 f and the provisional
attenuation bPif. Therefore, the modification of Eq. (9)assumes application of the formula expressed by
2pfdTi 2pfd 2pf ðdTi  dÞ
ﬃ di 00 f þ hi 00 ;cPi ð f Þ ci ð f Þ
cW

(11)

where di00 and hi00 are real constants. When the assumption
given by Eq. (11) is valid, the application of a time shift anda phase rotation to the conventional model can produce
appropriate fast and slow waves. The proposed method optimizes four parameters; the signal amplitude parameter Ai0 ,the temporal specimen thickness dTi, the time shift parameter
di0 , and the phase rotation hi0 .
Figure 2 shows the waveforms of the fast and slow
waves that were calculated using the modified transfer function in the experimental study, where the specimen thickness
varied from 5 to 25 mm and the phase rotation parameter
was 0 . The envelope amplitude of each wave was normalized, and the rise time of each wave was set to t ¼ 0 s after atime shift process. In this study, we defined the rise time ofeach wave as the time when the normalized amplitude firstexceeds 0.01, where specimen thicknesses are 5 and 20 mm.
The rise times for other thicknesses were determined byinterpolation and extrapolation. As shown in Fig. 2, the riseTaki et al.: Fast characterization of bone wave

1685

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21

FIG. 2. Waveforms of (a) fast and (b)slow waves calculated using the modified transfer function in the experimental study. The specimen thickness dranges from 5 to 25 mm. The phaserotation parameter hi is 0 .

times of the fast and slow waves were determined stably.
We can use another process to determine the rise time ofeach wave. The phase rotation parameter hi changes the
waveforms of the fast and slow waves, as shown in Fig. 3.
B. Deconvolution method based on the FDI imaging
methodThe proposed deconvolution method consists of twoparts; the first is the initial estimation of the fast and slow
waves using the FDI imaging method, and the second is theconclusive estimation using the least squares method in thetime domain. Figure 4 shows a flowchart for the proposed
deconvolution method.
1. Initial estimation using the FDI imaging method

Typically, the amplitude of the fast wave is lower thanthat of the slow wave. Because the fast wave is received earlier than the slow wave, the proposed method selects a time
window for estimation of a fast wave, where the response of
the slow wave is negligible within that time window. To
select a time threshold that separates the fast and slow
waves, a previous work used a local minimum in the envelope of the received signal in the time domain (Wear, 2014).
However, it is difficult to use this strategy under the strict
condition that the peak amplitude position of the fast wave islocated close to the fast wave in the time domain. The first
purpose of using the FDI imaging method is for accurate
selection of the rise time of the slow wave, TS. The second
purpose is to prepare candidate waves for the fast and slow
waves.

The FDI imaging method is a high-range-resolution
imaging method that uses an adaptive beamforming technique (Taki et al., 2012). The method suppresses the contributions of undesired echoes under the constraint that areference signal must give a constant response. Therefore,
we construct several candidate waves for the fast and slow
waves using the modified transfer function that was
described in Sec. II A, where each candidate is used as thereference signal in the FDI imaging method. The estimatedintensity calculated by the FDI imaging method is normalized with respect to the amplitude of the employed reference
signal. To eliminate the effect of this normalization on theestimated intensity, we set the peak amplitude of each envelope to 1.
The FDI imaging method was designed based on thecondition that a narrowband signal is used, and it assumes
that all frequency components of a single wave have thesame intensity and phase. During the application of ultrasound using a broadband signal, the intensity and the phaseof a frequency component of a single wave vary with the frequency. To ensure that the intensity and phase of the frequency components are uniform, each frequency componentof a received signal is normalized with respect to the reference signal (Taki et al., 2012),

XHl ¼ Xl XRl=ðjXRl j2 þ gÞ;(12)

where Xl, XHl, and XRl are the lth frequency components of
the received signal, of the received signal after normalization, and of the reference signal, respectively, the superscript
asterisk denotes the complex conjugate, and g is a constantterm used for stabilization. In this study, we set the value ofFIG. 3. Waveforms of (a) fast and (b)slow waves calculated using the modified transfer function in the experimental study. The specimen thickness d is
9 mm. The phase rotation parameters
hi are 0 , 45 , 90 , and 135 .

1686

J. Acoust. Soc. Am., Vol. 137, No. 4, April 2015

Taki et al.: Fast characterization of bone wave

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21

does not work when the waves are correlated with each other.
Therefore, we must suppress the correlation between the fastand slow waves. The employment of the frequency averaging
technique suppresses the correlation between the two waves(Taki et al., 2012),
RA ¼

M
1X
Rm ;
M m¼1

(16)

where RA is the covariance matrix after frequency averaging, M is the number of sub-matrices used for frequency
averaging, and Rm is a sub-matrix of R where the (i0 , j0 ) element of Rm corresponds to the (i0 þ m  1, j0 þ m  1) element of R.
The FDI imaging method with the Capon method suppresses any undesired echoes under the constraint of a constant response for the reference signal at a measurement
depth (Capon, 1969). This problem is expressed by the following formulas:min

WT RA W

C ¼ ½ ejk1 r

ejk2 rsubjecttoCT W ¼ 1 ;T   ejkL r  ;

(17)
(18)

where L is the size of RA. This problem can be solved by theapplication of Lagrange multiplier methods without the iteration process, as emphasized in this study (Taki et al., 2012).
When the optimum weighting vector is used, the estimatedintensity can be expressed using the following formulas:
PCap ðrÞ ¼

FIG. 4. Flowchart of the proposed deconvolution method.

g to be 40 dB lower than the average intensity of the reference signal, XRl.
The phase of each frequency component of a received
signal XHl depends on the product of its frequency and thepath length. The output intensity of the FDI imaging method
is therefore formulated as follows:PðrÞ ¼ jXTH W j2 ¼ WT RW;(13)

R ¼ XH XT
H ;

(14)

where r/2 is the measurement depth, XH is the received signal vector that consists of XHl, the superscript T denotes thetranspose, and W is a weighting vector used for phase compensation. The simplest setting is the employment of a fixedweighting vector given by
W ¼ ½ ejk1 r

ejk2 rT ;(15)

where kl is the lth wavenumber of the frequency components.
The proposed FDI imaging method uses the Capon
method, which is an adaptive beamforming technique.
However, the FDI imaging method using the Capon method
J. Acoust. Soc. Am., Vol. 137, No. 4, April 20151
C ðRA þ g0 EÞ1 C
T

;

(19)

where g0 E is a diagonal loading used for stabilization. In thisstudy, we set the value of g0 to be 40 dB lower than the average of the diagonal elements of RA. Equation (19) indicatesthat only a single calculation of the inverse matrix is
required to acquire the estimated intensity profile for all
depths without calculation of the optimum weighting vector.
The FDI imaging method is therefore suitable for real-time
signal processing applications.
When the waveform of a candidate is unlike that of theslow wave, the output intensity estimated by the FDI imaging method generally decreases. In other words, the output
intensity will be at a maximum when the candidate waveform is closest to the waveform of the slow wave. We therefore calculate the intensity profiles of all candidates for theslow waves, and select the candidate where the peak intensity in its profile is the highest among all candidates. We use
the depth of the peak intensity of the selected candidate forthe rise time of the slow wave. In the least mean squares process, we extracted the candidates with peak intensities thatwere higher than 1/100 of the maximum peak intensity. In
the FDI imaging process, we prepared 11 candidates foreach of the fast and slow waves, where the thickness dranged from 5 to 25 mm and the sampling interval is 2 mm.
In the fast wave case, the rise time and the candidate extraction for the least squares method were determined by thesame process that was used for the slow wave, with theTaki et al.: Fast characterization of bone wave

1687

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21

FIG. 5. Slow wave intensity estimatedusing the FDI imaging method when
the specimen thickness is (a) 9 mm and(b) 12 mm in the experimental study.
The vertical broken line corresponds to
the rise time of the slow wave as estimated by the FDI imaging method.

exception that the rise time of the fast wave should be atleast 0.4 ls earlier than that of the slow wave.
Figure 5 shows the intensity profile of the slow wavethat was estimated using the FDI imaging method in the experimental study. Wear used the local minimum of the envelope to determine the border between the fast and slow
waves (Wear, 2014). However, in many cases, no local minimum exists for the envelope, as shown in Fig. 5. In contrast,
the estimated intensity profile of the slow wave that was calculated by the FDI imaging method shows a sharp peak at
the rise time of the slow wave, TS. Therefore, the FDI imaging method is suitable for determination of the border
between the fast and slow waves.
2. Conclusive estimation using the least squares
methodThe proposed decomposition method uses the leastsquares method in the time domain to characterize the fastand slow waves. In the characterization process, the proposed method optimizes the signal amplitude parameter Ai0 ,the temporal specimen thickness dTi, the rise time
TRi ¼ di0 /2p, and the phase rotation hi0 , as expressed in Eq.
(9). When TRi ¼ 0, i.e., di0 ¼ 0, a candidate for the fast andslow waves in the time domain is expressed using the following formulas:si ðdTi ; Ai 0 ; hi 0 ; tÞ ¼ RefF1 ½SI ðf ÞHi 0 ðf Þdi 0 ¼0 g
ﬃ Ai 0 cos hi 0 Refs0i ðtÞgprovisional constant values for bPi and cPi, and the waveformof si0 (t) varies with the single variable dTi. Because eachvalue of Ai0 cos hi0 and Ai0 sin hi0 is a real number, estimationof Ai0 and hi0 requires a single linear least squares process on
two variables for each variable set of dTi and TRi. Therefore,
the number of linear least squares processes required is the
product of the number of specimen thickness candidates dTi
and the number of rise time candidates TRi.
In the least mean squares process, we prepared 111 candidates for each of the fast and slow waves, where d ranged
from 4 to 26 mm and the sampling interval was 0.2 mm.
Because candidate extraction using the FDI imaging method
reduced the number of candidates, the number of candidates
that was actually used in the least mean squares process
was 111 or less. The rise time candidates ranged from TiFDI
 0.8 ls to TiFDI þ 0.8 ls, where TiFDI is the rise time calculated using the FDI imaging method, and T2FDI ¼ TS.
For characterization of the fast and slow waves, weselected the fitting region TA t TP, as shown in Fig. 6. In
this study, TA ¼ TE  4 ls and TP ¼ TE þ 0.8 ls, where TE is
the peak time of the envelope. These settings were usedbecause the t < TA region contains almost no signal, and in
the region where t > TP, the contributions of multiple reflections in the cancellous bone might not be negligible.
In the least squares process, the fast and slow waves are
estimated alternately and recurrently. The objective functions used for fitting of the fast and slow waves are expressedby the following formula:(20)
þ Ai 0 sin hi 0 Refjs0i ðtÞg;




2pfdTi 2pfdTi
s0i ðtÞ ¼ F1 SI ð f Þexp bP i fdTi þj;
cPi ð f Þ
cW
(21)
where F1 denotes an inverse Fourier transform. Estimationof the fast and slow waves uses minimization of the sum ofsquares D on four variables,DðdTi ; Ai 0 ; hi 0 ; TRi Þ ¼

T2
X

jso ðtÞ  si ðdTi ; Ai 0 ; hi 0 ; t  TRi Þj2 ;

t¼T1

(22)
where so(t) is the objective function used for fitting, andT1 t T2 is a fitting region. The proposed method uses
1688

J. Acoust. Soc. Am., Vol. 137, No. 4, April 2015

FIG. 6. Fitting region used in the least squares method. TE is the peak timeof the envelope, and TS is the rise time of the slow wave calculated using
the FDI imaging method. TA ¼ TE  4 ls, and TP ¼ TE þ 0.8 ls.
Taki et al.: Fast characterization of bone wave

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21


so ðtÞ ¼

sR ðtÞ  sE2 ðtÞ ði ¼ 1Þ;sR ðtÞ  sE1 ðtÞ ði ¼ 2Þ;(23)

where sR(t) is the received signal, and sE1(t) and sE2(t) are
the fast and slow waves that were estimated in the previous
least mean squares process, respectively. At the beginning of
the estimation process, sE1(t) ¼ sE2(t) ¼ 0.
First, the slow wave is estimated using the least meansquares method, where the fitting region is set as TS t TP,
i.e., T1 ¼ TS and T2 ¼ TP in Eq. (22). To estimate the fast
wave, the method is applied to a region where the contribution of the slow wave is negligible. The fitting region for the
fast wave is therefore set as TA t TS. After estimation of
the fast wave, the slow wave is estimated for a fitting regionof TA t TP. The initial estimation of the slow wave uses
the posterior fitting region of TS t TP, because the fast
wave has not been eliminated from the objective function forslow wave fitting. The number of iteration steps at the fittingof the fast and slow waves ranged from 1 to 5.

decomposition method. Figure 8 shows the estimated fastand slow waves in the experimental study. In both cases, thesummation of the estimated fast and slow waves strongly
agreed with the received signal. Figure 9 shows the residual
intensity normalized with respect to the intensity of thereceived signal over the whole fitting region of TA t TP,
which has a width of 4.8 ls. The proposed method succeeded
to characterize of the fast and slow waves accurately, where
the residual intensity normalized with respect to the received
signal intensity was less than 20 dB in the experimental
study. In the simulation study, the residual intensity was lessthan 20 dB when the specimen thickness was from 3 to
8 mm. In both the simulation and experimental studies, thenumber of iteration steps had only a slight effect on the
results. These results demonstrate the high performance androbustness of the proposed method in the estimation of fastand slow waves. The calculation time required to estimate

III. RESULTS

Figure 7 shows the fast and slow waves in the simulation study that were estimated using the proposed

FIG. 7. Fast and slow waves estimated by the proposed decomposition
method for specimen thicknesses of (a) 4 mm, (b) 6 mm, and (c) 8 mm in thesimulation study. Each constructed signal represents the summation of theestimated fast and slow waves.
J. Acoust. Soc. Am., Vol. 137, No. 4, April 2015

FIG. 8. Fast and slow waves estimated by the proposed decomposition
method for specimen thicknesses of (a) 6 mm, (b) 9 mm, (c) 12 mm, and (d)15 mm in the experimental study. Each constructed signal represents thesummation of the estimated fast and slow waves.
Taki et al.: Fast characterization of bone wave

1689

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21

FIG. 9. Residual intensity in (a) thesimulation study and (b) the experimental study when normalized withrespect to the intensity of the received
signal over the entire fitting region,
which has a width of 4.8 ls. The number of iteration steps at the fitting offast and slow waves is from 1 to 5.

the fast and slow waves for a single received signal was 5 susing a laptop personal computer (PC) with a single central
processing unit (CPU), where the number of iteration steps
was 3.
IV. DISCUSSION

Figure 9 shows that a single estimation process was sufficient for accurate characterization of the fast and slow
waves, and that the number of iteration steps had only asmall effect on the characterization performance. The highperformance and robustness of the proposed method may
originate from the accurate estimation of the rise time of theslow wave, TS, using the FDI imaging method. Because thefitting region for the fast wave is t < TS, the accurate estimation of TS provides an appropriate fitting region for the fast
wave that is wide and is little affected by the estimated slow
wave. Therefore, the number of iteration steps may have asmall effect on the estimation of the fast wave using the proposed method, as indicated by the results shown in Fig. 9.
In most cases, the least mean squares process used all
111 candidates, especially for the slow wave fitting process.
In the simulation study, the average numbers of candidates
used by the process were 72.5 and 111 for the fast and slow
waves, respectively. In the experimental study, 89.2 and102.1 candidates were used for the fast and slow waves,respectively. The proposed method used the candidates withpeak intensities that were higher than 1/100 of the maximum
peak intensity, and a more severe constraint will suppress
the computational load, although it would be at the cost ofsome robustness.

The proposed decomposition method estimated thespecimen thickness, the signal amplitude and the phase rotation for each wave. We calculated the signal amplitudes Ai0
and the slopes of attenuation (broadband ultrasonic attenuation parameters) bi0 from the true specimen thickness andusing the parameters estimated by the proposed method, as
shown in Fig. 10. In the experimental study, the signal amplitude of the fast wave was higher than 1 under all specimen
thickness conditions. This is inconsistent with Eq. (5),because the signal amplitude should be less than 1. This
result indicates that the linear-with-frequency attenuation
model may be inappropriate for the fast wave transfer modelat least. The frequency components of the fast wave thatwere estimated by the proposed method may agree withthose of the true fast wave in the high signal-to-noise ratio
range, as shown in Fig. 11. Because the true signal amplitude
should be less than 1, the attenuation coefficient of the fast
wave a1(f) may be written using the following formula:a1 ðf Þ ¼ b1 f y

ðy > 1Þ:(24)

The use of an attenuation model described by a powerlaw, i.e., the model given by Eq. (24), should improve theestimation accuracy for characterization of the fast and slow
waves; however, this also causes a large increase in the computational load. Because the residual intensity normalizedwith respect to the received signal was less than 20 dB, thelinear-with-frequency attenuation model may be sufficientfor the proposed decomposition method.
Figure 12 shows the phase rotation values hi for the fastand slow waves that were estimated by the proposed method

FIG. 10. (a) Signal amplitudes and (b)slopes of attenuation (broadband ultrasonic attenuation parameters) for the
fast and slow waves in the experimental study estimated using the proposed
decomposition method with a linearwith-frequency attenuation model.

1690

J. Acoust. Soc. Am., Vol. 137, No. 4, April 2015

Taki et al.: Fast characterization of bone wave

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21

V. CONCLUSION

FIG. 11. Schema of the signal amplitude of a frequency component, where
the longitudinal axis uses a logarithmic scale. ATrue and AEst are the true signal amplitude and the estimated amplitude, respectively, and were obtainedusing the proposed method with the linear-with-frequency attenuation
model. The gray region denotes the frequency components with high signalto-noise ratios.

in the experimental study. The phase rotation of the slow
wave seemed to be constant. In contrast, the phase rotationof the fast wave varied, and its differential seemed to berelated to the differential of the specimen thickness. This
result can be explained by introducing the hypothesis thatwave transfer in a cancellous bone causes a physical phenomenon, e.g., the uneven wavefront, and the phase rotation
originates from this phenomenon. Because the slow wave
mainly transfers in the bone marrow, which is a liquid material, this physical phenomenon may not occur and the phaserotation may in fact be independent of the specimen thickness. The constant phase rotation may originate at the
moments of transmission and reception of the ultrasoundwave, where it might be caused by the waveform changethat originated from the uneven wavefront. In contrast, thefast wave mainly transfers in the bone structure. Therefore,
the physical phenomenon may occur at the point of transfer
in the bone, and the cumulative effect of this phenomenon
may cause the phase rotation that depends on the specimen
thickness. This consideration may validate the use of thephase rotation parameter in the transfer function.

FIG. 12. Phase rotation values hi for the fast and slow waves estimated using
the proposed method in the experimental study.
J. Acoust. Soc. Am., Vol. 137, No. 4, April 2015

In this study, we proposed a novel fast decomposition
method based on the FDI imaging method, where the FDI
imaging method enables accurate determination of the border between the fast and slow waves. In addition, we proposed a modified wave transfer function using a phaserotation parameter that may compensate for the waveformchange caused by the uneven wavefront. The proposed
method provided accurate characterization of the fast andslow waves, where the residual intensity normalized withrespect to the received signal intensity was found to be lessthan 20 dB in an experimental study using bone specimens
with thicknesses ranging from 6 to 15 mm. In the simulation
study, the residual intensity was less than 20 dB when thespecimen thickness was from 3 to 8 mm. The characterization of the fast and slow waves for a single received signal
takes only 5 s using a laptop PC with a single CPU. We
believe that this method has great potential to provide good
indicators of osteoporosis accurately and rapidly.
ACKNOWLEDGMENTS

This work was partly supported by the Innovative
Techno-Hub for Integrated Medical Bio-imaging Project of
the Special Coordination Funds for Promoting Science andTechnology from the Ministry of Education, Culture, Sports,
Science and Technology (MEXT), Japan, and by MEXT/JSPS KAKENHI through Grant Nos. 25870345 and
25871038.
Anderson, C. C., Bauer, A. Q., Holland, M. R., Pakula, M., Laugier, P.,
Bretthorst, G. L., and Miller, J. G. (2010). “Inverse problems in cancellous
bone: Estimation of the ultrasonic properties of fast and slow waves using
Bayesian probability theory,” J. Acoust. Soc. Am. 128, 2940–2948.
Anderson, C. C., Marutyan, K. R., Holland, M. R., Wear, K. A., and Miller,
J. G. (2008). “Interference between wave modes may contribute to theapparent negative dispersion observed in cancellous bone,” J. Acoust. Soc.
Am. 124, 1781–1789.
Barkmann, R., Laugier, P., Moser, U., Dencks, S., Klausner, M., Padilla, F.,Ha€ıat, G., and Gl€
uer, C. C. (2008a). “A device for in vivo measurementsof quantitative ultrasound variables at the human proximal femur,” IEEE
Trans. Ultrason. Ferroelectr. Freq. Control 55, 1197–1204.
Barkmann, R., Laugier, P., Moser, U., Dencks, S., Klausner, M., Padilla, F.,Ha€ıat, G., and Gl€
uer, C. C. (2008b). “In vivo measurements of ultrasoundtransmission through the human proximal femur,” Ultrasound Med. Biol.
34, 1186–1190.
Blake, G. M., and Fogelman, I. (2007). “Role of dual-energy X-ray absorptiometry in the diagnosis and treatment of osteoporosis,” J. Clin.
Densitometry 10, 102–110.
Bouxsein, M. L., Courtney, A. C., and Haynes, W. C. (1995). “Ultrasoundand densitometry of the calcaneus correlate with the failure loads of
cadaveric femurs,” Calcif. Tissue Int. 56, 99–103.
Capon, J. (1969). “High resolution frequency-wavenumber spectrum analysis,”
Proc. IEEE 57, 1408–1418.
Cardoso, L., Teboul, F., Sedel, L., Oddou, C., and Meunier, A. (2003). “Invitro acoustic waves propagation in human and bovine cancellous bone,”
J. Bone Miner. Res. 18, 1803–1812.
Dencks, S., and Schmitz, G. (2013). “Estimation of multipath transmission
parameters for quantitative ultrasound measurements of bone,” IEEE
Trans. Ultrason. Ferroelectr. Freq. Control 60, 1884–1895.
Fellah, Z. E., Chapelon, J. Y., Berger, S., Lauriks, W., and Depollier, C.
(2004). “Ultrasonic wave propagation in human cancellous bone:
Application of Biot’s theory,” J. Acoust. Soc. Am. 116, 61–73.
Genant, H. K., and Jiang, Y. (2006). “Advanced imaging assessment ofbone quality,” Ann. N.Y. Acad. Sci. 1068, 410–428.
Taki et al.: Fast characterization of bone wave

1691

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21

Ha€ıat, G., Padilla, F., Peyrin, F., and Laugier, P. (2007). “Variation of ultrasonic parameters with microstructure and material properties of trabecular bone: A 3D model simulation,” J. Bone Miner. Res. 22,
665–674.
Ha€ıat, G., Padilla, F., Svrcekova, M., Chevalier, Y., Pahr, D., Peyrin, F.,
Laugier, P., and Zysset, P. (2009). “Relationship between ultrasonic parameters and apparent trabecular bone elastic modulus: A numericalapproach,” J. Biomech. 42, 2033–2039.
Han, S., Medige, J., Davis, J., Fishkin, Z., Mihalko, W., and Ziv, I. (1997).
“Ultrasound velocity and broadband attenuation as predictors of loadbearing capacities of human calcanei,” Calcif. Tissue Int. 60, 21–25.
Hoffman, J. J., Nelson, A. M., Holland, M. R., and Miller, J. G. (2012).
“Cancellous bone fast and slow waves obtained with Bayesian probability
theory correlate with porosity from computed tomography,” J. Acoust.
Soc. Am. 132, 1830–1837.
Hosokawa, A. (2009). “Numerical analysis of variability in ultrasound propagation properties induced by trabecular microstructure in cancellous
bone,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 56, 738–747.
Hosokawa, A., and Otani, T. (1997). “Ultrasonic wave propagation in bovine cancellous bone,” J. Acoust. Soc. Am. 101, 558–562.
Hosokawa, A., and Otani, T. (1998). “Acoustic anisotropy in bovine cancellous bone,” J. Acoust. Soc. Am. 103, 2718–2722.
Hughes, E. R., Leighton, T. G., Petley, G. W., and White, P. R. (1999).
“Ultrasonic propagation in cancellous bone: A new stratified model,”Ultrasound Med. Biol. 25, 811–821.
Hughes, E. R., Leighton, T. G., White, P. R., and Petley, G. W. (2007).
“Investigation of an anisotropic tortuosity in a Biot model of ultrasonicpropagation in cancellous bone,” J. Acoust. Soc. Am. 121, 568–574.
Kaczmarek, M., Kubik, J., and Pakula, M. (2002). “Short ultrasonic waves
in cancellous bone,” Ultrasonics 40, 95–100.
Langton, C. M., Palmer, S. B., and Porter, R. W. (1984). “The measurementof broadband ultrasonic attenuation in cancellous bone,” Eng. Med. 13,
89–91.
Laugier, P. (2008). “Instrumentation for in vivo ultrasonic characterizationof bone strength,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 55,
1179–1196.
Lochm€uller, E. M., Eckstein, F., Zeller, J. B., Steldinger, R., and Putz, R.
(1999). “Comparison of quantitative ultrasound in the human calcaneus
with mechanical failure loads of the hip and spine,” Ultrasound Obstet.
Gynecol. 14, 125–133.

1692

J. Acoust. Soc. Am., Vol. 137, No. 4, April 2015

Marutyan, K. R., Bretthorst, G. L., and Miller, J. G. (2007). “Bayesian estimation of the underlying bone properties from mixed fast and slow mode
ultrasonic signals,” J. Acoust. Soc. Am. 121, EL8–EL15.
Marutyan, K. R., Holland, M. R., and Miller, J. G. (2006). “Anomalous negative dispersion in bone can result from the interference of fast and slow
waves,” J. Acoust. Soc. Am. 120, EL55–EL61.
Mizuno, K., Matsukawa, M., Otani, T., Laugier, P., and Padilla, F. (2009).
“Propagation of two longitudinal waves in human cancellous bone: Anin vitro study,” J. Acoust. Soc. Am. 125, 3460–3466.
Mizuno, K., Matsukawa, M., Otani, T., Takada, M., Mano, I., andTsujimoto, T. (2008). “Effects of structural anisotropy of cancellous boneon speed of ultrasonic fast waves in the bovine femur,” IEEE Trans.
Ultrason. Ferroelectr. Freq. Control 55, 1480–1487.
Nagatani, Y., Mizuno, K., Saeki, T., Matsukawa, M., Sakaguchi, T., andHosoi, H. (2008). “Numerical and experimental study on the wave attenuation in bone—FDTD simulation of ultrasound propagation in cancellous
bone,” Ultrasonics 48, 607–612.
Nelson, A. M., Hoffman, J. J., Anderson, C. C., Holland, M. R., Nagatani,Y., Mizuno, K., Matsukawa, M., and Miller, J. G. (2011). “Determining
attenuation properties of interfering fast and slow ultrasonic waves in cancellous bone,” J. Acoust. Soc. Am. 130, 2233–2240.
Njeh, C. F., Hans, D., Fuerst, T., Gl€uer, C. C., and Genant, H. K. (1999).
Quantitative Ultrasound: Assessment of Osteoporosis and Bone Status(Martin Dunitz, London), pp. 1–432.
Njeh, C. F., Kuo, C. W., Langton, M. M., Atrah, H. I., and Boivin, C. M.
(1997). “Prediction of human femoral bone strength using ultrasound velocity and BMD: An in vitro study,” Osteoporosis Int. 7, 471–477.
Taki, H., Taki, K., Sakamoto, T., Yamakawa, M., Shiina, T., Kudo, M., andSato, T. (2012). “High range resolution ultrasonographic vascular imaging
using frequency domain interferometry with the Capon method,” IEEE
Trans. Med. Imaging 31, 417–429.
Waters, K. R., Mobley, J., and Miller, J. G. (2005). “Causality-imposed(Kramers-Kronig) relationships between attenuation and dispersion,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52, 822–833.
Wear, K. A. (2013). “Estimation of fast and slow wave properties in cancellous bone using Prony’s method and curve fitting,” J. Acoust. Soc. Am.
133, 2490–2501.
Wear, K. A. (2014). “Time-domain separation of interfering waves in cancellous bone using band limited deconvolution: Simulation and phantom
study,” J. Acoust. Soc. Am. 135, 2102–2112.

Taki et al.: Fast characterization of bone wave

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 130.34.32.223 On: Mon, 27 Apr 2015 22:33:21