Fast decomposition of two ultrasound longitudinal waves incancellous bone using a phase rotation parameter for bonequality assessment: Simulation study
Hirofumi Takia)
Graduate School of Biomedical Engineering, Tohoku University, Sendai 980-8575, 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

Hiroshi Kanai
Graduate School of Engineering, Tohoku University, Sendai 980-8579, Japan

Shin-Ichi Izumi
Graduate School of Biomedical Engineering, Tohoku University, Sendai 980-8575, Japan(Received 3 May 2017; revised 31 August 2017; accepted 6 October 2017; published online 24
October 2017)
Ultrasound signals that pass through cancellous bone may be considered to consist of two longitudinal waves, which are called fast and slow waves. Accurate decomposition of these fast and slow
waves is considered to be highly beneficial in determination of the characteristics of cancellous
bone. In the present study, a fast decomposition method using a wave transfer function with a phaserotation parameter was applied to received signals that have passed through bovine bone specimens
with various bone volume to total volume (BV/TV) ratios in a simulation study, where the elastic
finite-difference time-domain method is used and the ultrasound wave propagated parallel to thebone axes. The proposed method succeeded to decompose both fast and slow waves accurately; thenormalized residual intensity was less than 19.5 dB when the specimen thickness ranged from 4
to 7 mm and the BV/TV value ranged from 0.144 to 0.226. There was a strong relationship between
the phase rotation value and the BV/TV value. The ratio of the peak envelope amplitude of thedecomposed fast wave to that of the slow wave increased monotonically with increasing BV/TV
ratio, indicating the high performance of the proposed method in estimation of the BV/TV value in
C 2017 Acoustical Society of America. https://doi.org/10.1121/1.5008502
cancellous bone. V[KAW]Pages: 2322–2331

I. INTRODUCTION

Bone strength is known to depend on not only bonequantity (bone mineral density), but also “bone quality”
which would be determined by several factors: microstructure, elasticity, and metabolism of bone. The gold standard
method for diagnosis of osteoporosis is dual x-ray absorptiometry. However, because this method measures the bonemineral density only, its performance in the estimation of
osteoporosis is limited (Genant and Jiang, 2006; Blake andFogelman, 2007). Quantitative ultrasound (QUS) is anothermodality that is used to evaluate bone quality by estimating
the speed of sound and the broadband ultrasonic attenuation
characteristics of the bone (Langton et al., 1984; Laugier,2008; Barkmann et al., 2008; Bouxsein et al., 1995; Njeh
et al., 1997; Han et al., 1997; Lochm€uller et al., 1999; Ha€ıat
et al., 2009). It may be considered that ultrasound signals
passing through cancellous bone consist of two longitudinal
waves, which are called fast and slow waves (Hosokawa andOtani, 1997, 1998; Hosokawa, 2010; Mizuno et al., 2008;a)Electronic mail: hirofumi.taki.a1@tohoku.ac.jp2322

J. Acoust. Soc. Am. 142 (4), October 2017

Mizuno et al., 2009; Ha€ıat et al., 2008; Padilla and Laugier,2000). Because the fast wave mainly propagates through thebone trabeculae, the properties of the wave reflect the structure of the cancellous bone (Hosokawa and Otani, 1997,
Mizuno et al., 2008; Mizuno et al., 2010). In addition,
another previous work reported that the amplitude of theslow wave is closely related to the bone volume (Otani
2005). Therefore, analysis of these two waves should reveal
both the state of the cancellous structure of the bone and thebone volume, which are useful in the evaluation of bonequality (Cardoso et al., 2003; Ha€ıat et al., 2007; Hosokawa
2009; Laugier and Ha€ıat, 2011). However, the fast and slow
waves overlap with each other in many cases, thus obstructing accurate analysis of the properties of the two waves.
Accurate decomposition of the fast and slow waves
would be highly beneficial for determination of the characteristics of cancellous bone. One strategy that has been used
to decompose these waves is application of a Bayesian
method that uses a Markov chain Monte Carlo method withsimulated annealing (Marutyan et al., 2007; Anderson et al.,2010; Nelson et al., 2011; Hoffman et al., 2012; Groopman
et al., 2015). The Bayesian decomposition method has

0001-4966/2017/142(4)/2322/10/$30.00

C 2017 Acoustical Society of AmericaVdemonstrated high performance in accurately characterizing
the fast and slow waves, but at the cost of a high computational load. Two other approaches that have been used arefast decomposition algorithms: the modified least-square
Prony method and the space-alternating generalized
expectation-maximization algorithm (Dencks and Schmitz
2013; Wear 2013, 2014). The Bayesian and the modified
least-square Prony methods yielded comparable results forthe ultrasonic properties of fast and slow waves in cancellous
bone; however, additional improvement in decomposition
performance is desirable.
We recently reported a fast decomposition method thatwas based on a minimum variance (MV) range beamformer,
called the MV beamformer for range resolution improvement (Taki et al., 2015). The method introduces a modified
wave transfer function with a phase rotation parameter thatcan compensate for the arrival time distribution that iscaused by path-length variations within a bone specimen.
The fast decomposition method has demonstrated high accuracy in characterization of the two waves on a similar levelto that of a Bayesian method, thus indicating the validity ofuse of the modified wave transfer function with the phaserotation parameter.
In the present study, we investigated the relationshipbetween the two decomposed waves and the bone volumeusing the above method to evaluate the performance of thedecomposition method for bone quality assessment. A previous study reported that the waveform of the ultrasound signal changed not only with the bone sample thickness but
also with the transducer size (Fujita et al., 2013). Therefore,
we applied the fast decomposition method with the modified
wave transfer function to received signals that had passed
through bone specimens with various bone volume to total
volume (BV/TV) ratios and transducer sizes in a simulation
study using the finite-difference time-domain (FDTD)
method, in which the bone sample was fully immersed inwater. The aims of the study can be considered to be twofold. First, we investigated the validity of the phase rotation
parameter that was used by the modified wave transfer function. Second, we evaluated the performance of the proposed
fast decomposition method with the modified wave transfer
function in estimation of the porosity of cancellous bone.

II. METHODS

The fast decomposition method uses a modified transfer
function with a phase rotation parameter. In this section, we
explain the simulation settings used to construct the received
signals that pass through bone specimens with various BV/
TV values and transducer sizes. We then briefly describe themodified transfer function and the decomposition strategy
that was reported in our previous work (Taki et al., 2015).
A. Simulation settings

To derive the signals of samples with various BV/TV
values, the elastic FDTD method was used. The FDTD simulation software was originally programmed by our group(Nagatani et al., 2006). For the FDTD simulations, a parallelepipedic bovine cancellous bone sample was used
(Nagatani et al., 2008). Three-dimensional bone models
were created using x-ray micro-focus computed tomography(x-ray lCT) images (SMX-100CT, Shimadzu, Kyoto, Japan)of an actual bone sample. The spatial resolution of the CT
images was 46.0 lm. Because the CT images are grayscale
images, the simulation models were created by selecting a
specific binarization threshold that determines the trabecular
thickness but does not influence the bone structure (Nagatani
and Tachibana, 2014). In this case, seven thresholds of 145,
146, 147, 148, 149, 150, and 151 on a 256-scale gray image
were selected, which meant that the BV/TV ratios of theconstructed models ranged from 0.144 to 0.226. Figure 1
shows the bone structure under the seven BV/TV conditions.
The bone model dimensions were 15  15  8 mm3 and thesimulation field dimensions were 16  16  13 mm3. To
investigate the gradual change of the waveform as it propagated inside the cancellous bone, the model was gradually
reduced in thickness from 8 to 4 mm at intervals of 1 mm.
The porous part of the cancellous bone was filled with water.
The experimental observed values for wave speed and attenuation in the bovine cortical bone were used for the trabecular part of the cancellous bone (Sasso et al., 2007, 2008;
Yamato et al., 2008; Yamamoto et al., 2008; Nagatani andTachibana, 2014). The parameters that were used in the simulations are shown in Table I. As the initial signal at the surface of the plane source, an experimentally observed

FIG. 1. Bone specimen structure usedin the simulation study.

J. Acoust. Soc. Am. 142 (4), October 2017

Taki et al.

2323

TABLE I. Parameters used in the FDTD simulations.
Material

Water

Trabecular boneDensity [g/cm3]
Velocity [m/s]

1.0
1483
-

2.0
4060
2450

Longitudinal wave
Shear wave

Hi 0 ðf Þ ¼ FTi ðf ÞHi ðf Þ
2pfd 2pfd
00
¼ Ai exp bi fd þ jþ jdi f þ jhi ;ci ðf Þ cW
(4)

waveform of a single sinusoidal sound wave at a frequencyof 1 MHz that passed through the water was used. The wave
propagation direction was parallel to the direction of the predominant bone axis. The received signal was then calculated
at the plane surface that was placed opposite the transmitter,
where the size of the receiving transducer was same as thatof the transmitting transducer. In the present study, we usedfour different sizes of transducers, 40%, 60%, 80%, and fullsize of the simulation field size, that is, we used square transducers with side lengths of 6.4, 9.7, 12.9, and 16 mm.

B. Modified wave transfer functions for fast and slow
waves

The modified wave transfer functions for the fast andslow waves are based on the wave transfer function that wasmodeled for waves passing through cancellous bone in previous studies (Marutyan et al., 2006; Anderson et al., 2008).
The 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 a 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 have introduced the assumption thatthe wavefront of an ultrasound wave passing through a bone
specimen in water is the same as that passing through wateronly pass. However, during ultrasound propagation throughcancellous bone, there will be multiple paths with various
path lengths, causing an uneven wavefront and complicatedchanges in the waveform. Because the ultrasound waves thatpass through these multiple paths of various lengths are typically received by a flat transducer, the path-length-variation(PLV) should affect the transfer function. The modified
wave transfer function approximates this effect using the following formula:
SR 0 ðf Þ ¼ SI ðf Þ½FT1 ðf ÞH1 ðf Þ þ FT2 ðf ÞH2 ðf Þ;(2)

FTi ðf Þ ﬃ ATi exp ½ ci f þ jðdi f þ hi Þ;(3)where FT1(f) and FT2(f) are the frequency components forthe fast and slow waves, respectively, that originate from the
PLV effect, and ATi, ci, di, and hi are real constants. cif andjdif denote the attenuation and the time shift, respectively. hi
is a phase rotation parameter that is independent of f. In alinear-with-frequency attenuation model, the modified wave
2324

J. Acoust. Soc. Am. 142 (4), October 2017transfer function that accounts for the PLV effect isexpressed using the following formula (Taki et al., 2015):

Ai 0 ¼ Ai ATi ;(5)cbi 0 ¼ bi þ i ;d(6)where Ai is a signal amplitude parameter that is independentof f, bi is the attenuation coefficient, d is the bone specimen
thickness, ci(f) is the phase velocity, and cW is the soundvelocity in water. For each specimen with constant thicknessd, Ai,0 and bi0 are real constants, and jdif denotes the time
shift. When the phase velocities of the two waves satisfy the
Kramers-Kronig relations, they can be expressed using thefollowing formula (Waters et al., 2005; Bauer et al., 2007;
O’Donnell et al., 1978):
 
1
1
bi
f

¼  2 ln;(7)pci ðf Þ ci ðf0 Þ
f0
where f0 is a reference frequency.
C. Fast deconvolution method based on MV range
beamformer

The fast deconvolution method consists of two steps:
the first step estimates the rise times of the fast and slow
waves using the MV range beamformer, and the second step
decomposes these two waves using the least squares methodin the time domain.
1. Rise time estimation using the MV range
beamformer

During estimation of the fast wave, it is preferable toselect a time window in which the contribution of the slow
wave is negligible. To determine the appropriate time window for characterization of the fast wave, the decomposition
method estimates the slow wave rise time using the MV
range beamformer with frequency domain interferometry.
We prepared multiple candidate waves for use as the fastand slow waves using the modified wave transfer function.
The MV range beamformer uses each of the candidate waves
as a reference wave, and it then calculates the rise time andthe estimated intensity for each candidate wave separately.
The received signal is then normalized using a reference signal to adjust the intensities and phases of all frequency components uniformly (Taki et al., 2012; Taki et al., 2015):
XHl ¼ Xl XRl  =ðjXRl j2 þ gÞ;(8)where Xl, XHl, and XRl are the lth frequency components of thereceived signal, of the received signal after normalization, andof the reference signal, respectively; the superscript []* denotes
a complex conjugate, and g is a constant term that is used for
Taki et al.

stabilization. Because the optimum filter is the Wiener filter,the optimum value of eta should be equal to the noise intensity
for each frequency component. When the signal intensity isconstant and only white noise exists, we can estimate the noise
intensity of each frequency component using the varianceexpectation of the intensity of the component. However, in
medical ultrasound imaging, the signal intensity is commonly
variable and speckle is signal correlated noise. In the present
study, we followed our previous study (Taki et al., 2015) andintroduced the assumption that the expectation of the noise
intensity is 40 dB lower than the average intensity of the reference signal, that is, g was set to be 40 dB lower than the average intensity of the reference signal XRl.
The MV range beamformer minimizes the output intensity
under the constraint of a constant response for the reference
signal at a specific measurement depth (Taki et al., 2012). This
problem can be expressed using the following formulae:min

WT RA W

RA ¼

M
1X

Xm XTm ;M m¼1

subjectto CT W ¼ 1 ;

(9)
(10)


T
Xm ¼ XHm XHðmþ1Þ    XHðmþL1Þ ;
TC ¼ ejx1 t ejx2 t    ejxL t ;(11)where RA is the covariance matrix of the normalized frequency components after frequency averaging, W is aweighting vector, xl is the lth angular frequency, L is thesize of RA, and M is the number of sub-matrices that wereused for frequency averaging. The estimated intensity using
the optimum weighting vector can be expressed as
PMVR ðtÞ ¼

1
T

C ðRA þ g0 EÞ1 C

;

(12)

where g0 E is a diagonal loading that was used for stabilization. In the present study, the value of g0 was set to be 40 dB
lower than the average value of the diagonal elements of RA.
The size of M is set to be 40% of the total number of frequency components.
We calculate the intensity profiles PMVR(t) for all the
candidate reference waves for the slow waves, and thenselect the time at which the intensity is the highest among all
the candidates. The selected time is the tentative rise time ofthe slow wave TS. The decomposition procedure of the fast
wave using the least squares method uses the received signal
at which t < TS to eliminate the slow wave’s contribution.
We prepared 11 candidate reference waves for each of thefast and slow waves, where the thickness d in Eq. (4) ranged
from 5 to 25 mm with a sampling interval of 2 mm. We usedcandidates with peak intensities that were higher than 1/100of the maximum peak intensity for the fitting procedureusing the least squares method.
2. Decomposition using the least squares method

The proposed decomposition method uses the leastsquares method in the time domain (Taki et al., 2015). To
J. Acoust. Soc. Am. 142 (4), October 2017suppress the computational complexity of the process, we
modify the wave transfer function of Eq. (4) using propagation parameters that were reported in a previous study, as
described below. The wave transfer function used in thedecomposition process is given by
2pfdTi 2pfdTi
Hi 0 ðf Þ ﬃ Ai 0 exp bPi fdTi þ j
cPi ðf ÞcW
(13)þjdi 0 f þ jhi 0 ;
 
1
1
b i
f
;

¼  P2 lnp
cPi ðf Þ ci ðf0 Þ
f0

(14)

where dTi is the temporal specimen thickness, and bPif andcPi(f) are the provisional attenuation coefficient and phasevelocity that were reported in the previous study (Nelson
et al., 2011). The proposed method optimizes the following
four parameters for each wave: the signal amplitude parameter Ai0 , the temporal specimen thickness dTi, the time shiftparameter di0 , and the phase rotation hi0 . Because the proposed method characterizes each wave alternately, the fourparameters are all estimated simultaneously.
When the rise time TRi ¼ 0, candidates for the fast andslow waves in the time domain are expressed using the following formulae:nhiosi ðdTi ; Ai 0 ; hi 0 ; tÞ ¼ Re F1 SI ðf ÞHi 0 ðf ÞTRi ¼0
ﬃ Ai 0 cos hi 0 Refsi 0 ðtÞg
þ Ai 0 sin hi 0 Refjsi 0 ðtÞg;
2pfdTisi 0 ðtÞ ¼ F1 SI ðf Þexp bP i fdTi þ j
cPi ðf Þ
2pfdTiþ jdi 00 f ;

cW

(15)

(16)

where F–1 denotes an inverse Fourier transform. di00 is usedin the setting of TRi ¼ 0. We prepare si0 (t) for all candidate
waves in advance for use by the proposed method in thedecomposition process.

FIG. 2. Fitting region used in decomposition of fast and slow waves using
the least squares method. TE is the peak time of the envelope, and TS is therise time of the slow wave that was calculated using the frequency domaininterferometry (FDI) imaging method. TA ¼ TE  6 ls, and TP ¼ TE þ 0.8 ls.
Taki et al.

2325

The proposed method characterizes the fast and slow
waves alternately based on minimization of the sum ofsquares D of four variables:
T2
XDðdTi ;A0i ;h0i ;TRi Þ¼ jso ðtÞsi ðdTi ;A0i ;h0i ;tTRi Þj2 ; (17)t¼T1

where so(t) is the objective function used for the fitting,and T1  t  T2 is the fitting region. Because Ai0 coshi0 and
Ai0 sinhi0 are real numbers, the estimation of Ai0 and hi0 only
requires a single linear least squares process. Therefore,
the required number of linear least squares processes isgiven by the product of the number of candidate specimen

FIG. 3. Fast and slow waves that were estimated using the proposed decomposition method for specimen thicknesses of 4, 6, and 8 mm in the simulation study.
Each constructed signal represents the summation of the estimated fast and slow waves, where the bone volume to total volume (BV/TV) ratio ranged from
0.144 to 0.226 and the transducer size was 16 mm.
2326

J. Acoust. Soc. Am. 142 (4), October 2017

Taki et al.

thicknesses dTi and the number of candidate rise times TRi.
The objective function so(t) for the fast and slow waves isexpressed usingsR ðtÞ  sE2 ðtÞ ðfast waveÞ;(18)
so ðtÞ ¼
sR ðtÞ  sE1 ðtÞ ðslow waveÞ;where sR(t) is the received signal in the time domain, andsE1(t) and sE2(t) are the fast and slow waves that were estimated using the previous least mean squares process,respectively.
We prepared 111 candidates for each of the fast andslow waves, where d ranged from 4 to 26 mm using a sampling interval of 0.2 mm. The candidate rise times ranged
from TiMVR – 0.8 ls to TiMVR þ 0.8 ls for the fast and slow
waves, where TiMVR is the rise time that was calculated using
the MV range beamformer. Figure 2 shows the fitting region
that was used in the decomposition process using the leastsquares method. The fitting region in the present study
ranged from TA ¼ TE  6 ls to TP ¼ TE þ 0.8 ls, where TE is
the peak time of the envelope. We followed our previous
study to set these parameters (Taki et al., 2015). The t < TAregion contains almost no signal. In the region where t > TP,
the contributions of multiple reflections in the cancellous
bone might not be negligible. First, the proposed methodcharacterizes the slow wave using the received data within afitting region of TS  t  TP, because the slow wave rise timeshould be close to the rise time TS that was estimated using
the MV range beamformer. During the fast wave characterization, the proposed method uses the received data within afitting region of TA  t  TS to eliminate the contribution of
the slow wave. Following characterization of the fast wave,
the proposed method characterizes the slow wave using thereceived data from the whole fitting region, represented byTA  t  TP.

16 mm (Nagatani and Tachibara, 2014). The signals thatwere constructed by the summation of the decomposed fastand slow waves agreed strongly with the received signals.
Figure 4 shows the residual intensity when normalized withrespect to the received signal intensity over the fitting region.
The normalized residual intensity was less than 19.5 dB,
except under the conditions where the specimen thickness
was 8 mm and the BV/TV ratio was 0.211 or more. Figure 5
shows the fast and slow waves that were decomposed using
the proposed method, where the BV/TV ratio was 0.184 andthe transducer size ranged from 6.4 to 12.9 mm. Figure 6
shows the residual intensity when normalized with respect tothe received signal intensity over the fitting region, where
the normalized residual intensity was less than 22.7 dB.
These results demonstrate the excellent performance of the
proposed method in decomposition of the fast and slow
waves. The decomposition of two waves requires 5.5 s usinga Laptop PC with a single central processing unit (CPU).
Figure 7 shows the peak amplitudes of the envelopes ofthe fast and slow waves that were decomposed by the proposed method, where the BV/TV ratio ranged from 0.144 to
0.226 and the transducer size was 16 mm. The amplitude of
the fast wave envelope increased as the BV/TV value
increased; in contrast, that of the slow wave decreased under

III. RESULTS

Figure 3 shows the fast and slow waves that were
decomposed using the proposed method, where the BV/TV
ratio ranged from 0.144 to 0.226 and the transducer size was

FIG. 4. Residual intensity normalized with respect to the received signal
intensity over the fitting region for specimen thicknesses ranging from 4 to
8 mm, where the BV/TV ratio ranged from 0.144 to 0.226 and the transducer
size was 16 mm.
J. Acoust. Soc. Am. 142 (4), October 2017

FIG. 5. Fast and slow waves that were estimated by the proposed decomposition method for a specimen thickness of 4 mm in the simulation study.
Each residual signal was obtained by subtracting the estimated fast and slow
waves from the received signal, where the BV/TV value was 0.184 and thetransducer size ranged from 6.4 to 12.9 mm.
Taki et al.

2327

FIG. 6. Residual intensity when normalized with respect to the received signal intensity over the fitting region for specimen thicknesses ranging from 4
to 8 mm, where the BV/TV value was 0.184 and the transducer size ranged
from 6.4 to 16 mm.

the same conditions. Therefore, the ratio of the peak envelope amplitude of the fast wave to that of the slow wave
increases monotonically, as shown in Fig. 8. Because theincrease in the BV/TV value indicates thicker trabeculae,
this result is consistent with the consensus that the fast waves
mainly propagate in the trabeculae. Figure 9 shows the ratioof the peak envelope amplitude of the fast wave to that of
the slow wave, where the BV/TV ratio was 0.184 and thetransducer size ranged from 6.4 to 16 mm. The amplitude
ratio decreased as the transducer size increased from 6.4 to
16 mm; however, the observed variation was much smaller
than that caused by a change in the BV/TV ratio from 0.144
to 0.226, and the amplitude ratio of BV/TV ¼ 0.184 whenusing 6.4 mm transducers was close to that for BV/
TV ¼ 0.198 when using 16 mm transducers.
Figure 10 shows the phase rotation parameters for the
fast and slow waves that were estimated by the proposed
method, where the BV/TV ratio ranged from 0.144 to 0.226
and the transducer size was 16 mm. In both the fast and slow
waves, the phase rotation variations with respect to BV/TV
increased in magnitude with increasing bone specimen thickness. The phase rotation value of the estimated fast wave

FIG. 7. Peak amplitudes of the envelopes of the fast and slow waves thatwere decomposed by the proposed method for specimen thicknesses ranging
from 4 to 8 mm, where the BV/TV value ranged from 0.144 to 0.226 and thetransducer size was 16 mm.
2328

J. Acoust. Soc. Am. 142 (4), October 2017

FIG. 8. Ratio of the peak envelope amplitude of the fast wave to that of theslow wave for specimen thicknesses ranging from 4 to 8 mm, where the BV/
TV value ranged from 0.144 to 0.226 and the transducer size was 16 mm.

decreased as the BV/TV value increased. In contrast, the corresponding value of the slow wave increased. The large
phase rotation values and small residual intensity shown in
Fig. 4 illustrate the necessity of using the phase rotation
parameter in the wave transfer function.
Figure 11 shows the phase rotation parameters for the
fast and slow waves that were estimated by the proposed
method, where the BV/TV ratio was 0.184 and the transducer size ranged from 6.4 to 16 mm. The phase rotation variation caused by the change in the transducer size was muchsmaller than that which was caused by changes in the BV/
TV value. This result indicates the relationship that exists
between the BV/TV ratio and the phase rotations of the two
waves.
IV. DISCUSSION

This study employed an image processing algorithm tomodify bone volume fraction by using the same x-ray lCT
images of a sample (Ha€ıat et al., 2007; Nagatani et al.,2017). The binarizing threshold applied to the original x-ray
lCT images was varied so that the bone volume fraction of
the model was varied. We should notice that, by changing

FIG. 9. Ratio of the peak envelope amplitude of the fast wave to that of theslow wave for specimen thicknesses ranging from 4 to 8 mm, where the BV/
TV value was 0.184 and the transducer size ranged from 6.4 to 16 mm.
Taki et al.

the threshold, the connectivity of the trabeculae might beslightly changed and the direction of the thickness change
might not be exactly omnidirectional. This phenomenon canoccur in the in vivo situation, e.g., osteoporosis. Future work
should include the investigation of the limitation of thisimage processing algorithm.
The ratio of the peak envelope amplitude of the fast
wave to that of the slow wave increased monotonically withincreasing BV/TV, as shown in Fig. 8. The amplitude ratio
increased as the transducer size decreased; however, a 40%
reduction in the transducer size caused an increase in amplitude ratio similar to that caused by a 10% increase in BV/
TV. Because the appearance of the two-wave phenomenon
was reported in the case of waves passing through cancellous
bone with cortical bone (Nagatani et al., 2014), these results
indicate that accurate decomposition of the fast and slow
waves will enable precise estimation of the value of the BV/
TV ratio in cancellous bone, and that the proposed decomposition method is highly promising as a novel technique forosteoporosis diagnosis using QUS.
Figure 10 shows the strong relationship that exists
between the phase rotation value and the value of the BV/
TV ratio, where variation of the phase rotation increases
with increasing specimen thickness. The increase in the BV/
TV value indicates thicker bone trabeculae in the cancellous

FIG. 11. Phase rotation parameters for (a) fast and (b) slow waves for specimen thicknesses ranging from 4 to 8 mm, where the BV/TV value was 0.184
and the transducer size ranged from 6.4 to 16 mm.

FIG. 10. Phase rotation parameters for (a) fast and (b) slow waves for specimen thicknesses ranging from 4 to 8 mm, where the BV/TV value ranged
from 0.144 to 0.226 and the transducer size was 16 mm.
J. Acoust. Soc. Am. 142 (4), October 2017

bone, and the variation in the path length increases withincreasing specimen thickness. Because both the increase ininteraction through the surfaces of the bone trabeculae andthe increase in the specimen thickness should lead to variation in the path length distribution, this result is consistent
with the results of the method in which the phase rotation
originates from the path length distribution. The phase rotation of the fast wave decreased as the value of BV/TV
increased. In contrast, the phase rotation of the slow wave
increased under the same conditions. This phenomenon may
also have originated from the path-length variations.
Both the amplitude ratio and the phase rotation were
almost independent of the specimen thickness when the BV/
TV ratio was approximately 0.17, as shown in Figs. 8 and10. This phenomenon should be investigated in future work;
however, it may simply be caused by the balance of theultrasound signals that pass through the trabecular surface.
When the contribution of the ultrasound that goes from thetrabecular bone to the bone marrow is equivalent to thatfrom the bone marrow to the trabecular bone at each frequency, the contributions of the ultrasound signals that passthrough the trabecular surface may be negligible, thus suppressing the effects of the specimen thickness on the signal
waveform.
Taki et al.

2329

V. CONCLUSION

In the present study, we apply a fast decomposition
method that uses the wave transfer function with a phaserotation parameter to received signals that have passed
through bone specimens with various BV/TV values in asimulation study. The normalized residual intensity was lessthan 19.5 dB, except under the conditions that the specimen thickness was 8 mm and the BV/TV ratio was 0.211 or
more, thus demonstrating the high performance of the proposed method in decomposition of the fast and slow waves.
The strong relationship between the phase rotation value andthe value of BV/TV validated the use of the phase rotation
parameter in the wave transfer function to compensate forwaveform changes caused by the PLV effect. The ratio of
the peak envelope amplitude of the decomposed fast wave tothat of the slow wave increased monotonically with increasing BV/TV, thus demonstrating the high performance of the
proposed method in estimation of the value of BV/TV incancellous bone.
ACKNOWLEDGMENTS

This work was supported in part by the Life Science
Support Program of the Terumo Foundation for Life Scienceand Arts, by a research grant from the Suzuken Memorial
Foundation, and by KAKENHI (Grant No. 16K01431) fromthe Japan Society for the Promotion of Science (JSPS).
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. (2008). “A device for in vivo measurements ofquantitative ultrasound variables at the human proximal femur,” IEEE
Trans. Ultrason. Ferroelectr. Freq. Control 55, 1197–1204.
Bauer, A. Q., Marutyan, K. R., Holland, M. R., and Miller J. G. (2007). “Is
the Kramers-Kronig relationship between ultrasonic attenuation and dispersion maintained in the presence of apparent losses due to phase cancellation?,” J. Acoust. Soc. Am. 122, 222–228.
Blake, G. M., and Fogelman, I. (2007). “Role of dual-energy x-ray absorptiometry in the diagnosis and treatment of osteoporosis,” J. Clin. Densitom.
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.
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.
Fujita, F., Mizuno, K., and Matsukawa, M. (2013). “An experimental study
on the ultrasonic wave propagation in cancellous bone: Waveform changes
during propagation,” J. Acoust. Soc. Am. 134, 4775–4781.
Genant, H. K., and Jiang, Y. (2006). “Advanced imaging assessment ofbone quality,” Ann. N.Y. Acad. Sci. 1068, 410–428.
Groopman, A. M., Katz, J. I., Holland, M. R., Fujita, F., Matsukawa, M.,
Mizuno, K., Wear, K. A., and Miller, J. G. (2015). “Conventional,
Bayesian, and Modified Prony’s methods for characterizing fast and slow
waves in equine cancellous bone,” J. Acoust. Soc. Am. 138, 1594–1604.
2330

J. Acoust. Soc. Am. 142 (4), October 2017

Ha€ıat, G., Padilla, F., Peyrin, F., and Laugier, P. (2007). “Variation of ultrasonic parameters with microstructure and material properties of trabecularbone: A 3D model simulation,” J. Bone Miner. Res. 22, 665–674.
Ha€ıat, G., Padilla, F., Peyrin, F., and Laugier, P. (2008). “Fast wave ultrasonic propagation in trabecular bone: Numerical study of the influence ofporosity and structural anisotropy,” J. Acoust. Soc. Am. 123, 1694–1705.
Ha€ıat, G., Padilla, F., Svrcekova, M., Chevalier, Y., Pahr, D., Peyrin, F.,
Laugier, P., and Zysset, P. (2009). “Relationship between ultrasonicparameters 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. (2010). “Effect of porosity distribution in the propagation
direction on ultrasound waves through cancellous bone,” IEEE Trans.
Ultrason. Ferroelectr. Freq. Control 57, 1320–1328.
Hosokawa, A., and Otani, T. (1997). “Ultrasonic wave propagation inbovine 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.
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.
Laugier, P., and Haiat G., eds. (2011). Bone Quantitative Ultrasound(Springer, Dordrecht, Netherlands), pp. 1–6.
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.
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.
Mizuno, K., Somiya, H., Kubo, T., Matsukawa, M., Otani, T., andTsujimoto, T. (2010). “Influence of cancellous bone microstructure on two
ultrasonic wave propagations in bovine femur: An in vitro study,”
J. Acoust. Soc. Am. 128, 3181–3189.
Nagatani, Y., Guipieri, S., Nguyen, V. H., Chappard, C., Geiger, D., Naili,S., and Ha€ıat, G. (2017). “Three–dimensional simulation of quantitative
ultrasound in cancellous bone using the echographic response of a metallic
pin,” Ultrason Imag. 39, 295–312.
Nagatani, Y., Imaizumi, H., Fukuda, T., Matsukawa, M., Watanabe, Y., andOtani, T. (2006). “Applicability of finite-difference time-domain method
to simulation of wave propagation in cancellous bone,” Jpn. J. Appl. Phys.
45, 7186–7190.
Nagatani, Y., Mizuno, K., and Matsukawa, M. (2014). “Two-wave behavior
under various conditions of transition area from cancellous bone to cortical bone,” Ultrasonics 54, 1245–1250.
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.
Nagatani, Y., and Tachibana, R. O. (2014). “Multichannel instantaneous frequency analysis of ultrasound propagating in cancellous bone,” J. Acoust.
Soc. Am. 135, 1197–1206.
Taki et al.

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., 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.
O’Donnell, M., Jaynes, E., and Miller, J. G. (1978). “General relationships
between ultrasonic attenuation and dispersion,” J. Acoust. Soc. Am. 63,
1935–1937.
Otani, T. (2005). “Quantitative estimation of bone density and bone quality
using acoustic parameters of cancellous bone for fast and slow waves,”
Jpn. J. Appl. Phys. 44, 4578–4582.
Padilla, F., and Laugier, P. (2000). “Phase and group velocities of fast andslow compressional waves in trabecular bone,” J. Acoust. Soc. Am. 108,
1949–1952.
Sasso, M., Ha€ıat, G., Yamato, Y., Naili, S., and Matsukawa, M. (2007).
“Frequency dependence of ultrasonic attenuation in bovine cortical bone:
An in vitro study,” Ultrasound Med. Biol. 33, 1933–1942.
Sasso, M., Ha€ıat, G., Yamato, Y., Naili, S., and Matsukawa, M. (2008).
“Dependence of ultrasonic attenuation on bone mass and microstructure inbovine cortical bone,” J. Biomech. 41, 347–355.

J. Acoust. Soc. Am. 142 (4), October 2017

Taki, H., Nagatani, Y., Matsukawa, M., Mizuno, K., and Sato T. (2015). “Fastcharacterization of two ultrasound longitudinal waves in cancellous bone using
an adaptive beamforming technique,” J. Acoust. Soc. Am. 137, 1683–1692.
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. Imag. 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.
Yamamoto, K., Yaoi, Y., Yamato, Y., Yanagitan, T., Matsukawa, M., andYamazaki, K. (2008). “Ultrasonic wave properties in bone axis directionof bovine cortical bone,” Jpn. J. Appl. Phys. 47, 4096–4100.
Yamato, Y., Matsukawa, M., Yanagitani, T., Yamazaki, K., Mizukawa, H.,
and Nagano, A. (2008). “Correlation between hydroxyapatite crystalliteorientation and ultrasonic wave velocities in bovine cortical bone,” Calcif.
Tissue Int. 82, 162–169.

Taki et al.

2331