. ‘
u '"m‘n‘ﬁlﬂ‘h’ w z" ‘
a 4 C ‘ LIBRARY
M1 I Michg.. State
University
This is to certify that the
dissertation entitled
A VISCOELASTIC FINITE DIFFERENCE TIME DOMAIN
MODEL OF HUMAN THORAX TO DEVELOP AND VALIDATE
SOURCE LOCALIZATION ALGORITHMS
presented by
Sridhar Ramakrishnan
has been accepted towards fulﬁllment
of the requirements for the
Doctoral degree in Electrical Engineering
ML W
Major ProfessorIs Signature
limit 7, WWI
Q ﬁ
Date
MSU is an Aﬁ‘innative Action/Equal Opportunity Employer
— —.o.._.—r—uu.oonnonIllc
PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.
DATE DUE DATE DUE DATE DUE
Wu
5/08 K:lProlecc&Pres/ClRC/DateDue‘indd
A VISCOELASTIC FINITE DIFFERENCE TIME DOMAIN MODEL OF
HUMAN THORAX TO DEVELOP AND VALIDATE SOURCE
LOCALIZATION ALGORITHMS
By
Sridhar Ramakrishnan
A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Electrical Engineering
2009
ABSTRACT
A VISCOELASTIC FINITE DIFFERENCE TIME DOMAIN MODEL OF HUMAN
THORAX TO DEVELOP AND VALIDATE SOURCE LOCALIZATION
ALGORITHMS
By
Sridhar Ramakrishnan
Auscultation sounds offer a rich source of diagnostic information that could
potentially be used to detect a vast number of heart and lung pathologies noninvasively.
Recent times have seen efforts directed towards the possibility of using an array of
acoustic sensors for the localization of sounds, thereby leading to improved diagnosis.
However, a majority of the existing array processing and inversesource solution schemes
proposed to date rely on simplifying assumptions, and many of the factors that contribute
to violation of the assumptions, such as heterogeneity of the thoracic cavity, shear wave
contributions, extreme nearﬁeld conditions, to name a few, are not taken into account.
This work presents a test bed capable of simulating the sound distribution around the
thorax. Speciﬁcally, a ﬁnitedifference timedomain (FDTD) forward model of the
human thorax that solves the viscoelastic wave equations while taking the intrinsic
anatomy and the associated viscoelastic properties of the various tissues and structures
into account is presented. The test bed can be employed to develop and validate various
array signal processing algorithms and inverse problem techniques used for localizing
sound. Test results demonstrating the utility as well as effectiveness of the test bed are
presented.
ACKNOWLEDGMENTS
I wish to express my sincere gratitude to my advisor Dr. Satish Udpa and my
committee members Dr. Lalita Udpa, Dr. Joe Hauptman, and Dr. Ramakrishna
Mukkamala for their guidance and support in my PhD study. Their advice on technical
matters, encouragement, motivation and constructive criticism were greatly
acknowledged. I owe special thanks to Dr. Satish Udpa for having bolstered my
conﬁdence at every stage of my graduate student life and guiding me in making wise
career decisions.
Throughout the period of my graduate study, I was ﬁnancially supported by various
research projects undertaken by Dr. Satish Udpa, Dr. Lalita Udpa and Dr. Pradeep
Ramuhalli. This support, specially the one provided by Electric Power Research Institute,
is gratefully acknowledged. Being a member of the Nondestructive Evaluation
Laboratory, I feel gifted having had an opportunity to work with a team of dedicated
graduate students, interns, visiting scholars and professors over the past 6 years on
different research problems, which have not only widened my academic experience but
have also improved my interpersonal skills and instilled team ethics in me. I sincerely
thank all my laboratory members for thus having enriched my graduate student
experience. I also acknowledge Prof. Antonello Tambunino for his tutoring in the Inverse
Problems course and his valuable advice during my PhD.
I would like to thank the Department of Electrical and Computer Engineering, the
High Performance Computing Center and the College of Engineering at Michigan State
University for providing the best research resources and facilities. I thank the entire staff
iv
of Electrical and Computer Engineering, Office of the Dean of College of Engineering,
0188, DECS and Graduate School for their valuable timely assistance on various
occasions. I would like to specially acknowledge Linda Clifford for not only having
assisted me in all the paper work throughout the course of my graduate studies but also
for helping all of our laboratory members to forge familylike bonds amongst our group.
Finally, I wish to express my heartfelt gratitude to my parents, my sister and my
brotherin—law for their constant support and encouragement. Without their support this
dissertation would not have been possible.
TABLE OF CONTENTS
LIST OF TABLES ........................................................................................................... viii
LIST OF FIGURES ............................................................................................................. x
NOMENCLATURE .......................................................................................................... xv
CHAPTER 1 ......................................................................................................... 1
INTRODUCTION ............................................................................................................... 1
1.1 General ....................................................................................................................... 1
1.2 Dissertation layout ..................................................................................................... 3
CHAPTER 2 ......................................................................................................... 4
PRIOR WORK .................................................................................................................... 4
2.1 Introduction ............................................................................................................... 4
2.2 Acoustic Modeling Studies of the Thorax ...................... . .......................................... 5
2.3 Signal Processing Methods for Acoustic Diagnostics ............................................. 10
CHAPTER 3 ....................................................................................................... 12
VISCOELASTIC FDTD MODEL FOR HUMAN THORAX .............. 12
3.1 Introduction ............................................................................................................. 12
3.2 Governing Equations ............................................................................................... 12
3.2.1 Elastodynamic Equations ................................................................................... 13
3.2.2 Popular \flscoelastic Models .............................................................................. 15
3.2.3 Viscoelastic Wave Propagation Equation ........................................................ 16
3.3 Thorax Anatomy and Material Properties ............................................................... 19
3.4 Boundary Conditions ............................................................................................... 25
3.4.1 Concept of PML .................................................................................................... 25
3.4.2 Recursive Integration PML ................................................................................. 26
3.5 Source Model and Initial Conditions ....................................................................... 30
3.6 FiniteDifference TimeDomain Implementation ................................................... 33
3.6.1 Discretization scheme ......................................................................................... 33
3.6.2 Numerical Stability ............................................................................................... 38
3.7 Numerical Results .................................................................................................... 40
3.7.1 Validation of FDTD codes ................................................................................... 40
3.7.1.1 Comparison with analytical solution ................................................................ 40
3.7.1.2 Comparison with results published in literature ............................................... 46
3.7.2 Simulation results on human thorax .................................................................. 48
3.7.3 Wave Propagation Characteristics .................................................................... 57
3.8 Effects of Anatomical Variations ............................................................................ 60
3.8.1 Effects of ribs and shear elasticity of tissues ................................................... 60
3.8.2 Effect of size of thorax ......................................................................................... 62
CHAPTER 4 ........................................................................................... 64
STUDY OF SOUNDS UNDER PNEUMOTHORAX CONDITIONS ............................ 64
4.1 Introduction ............................................................................................................. 64
4.2 Simulation of Pneumothorax ................................................................................... 66
4.3 Results and Conclusion ........................................................................................... 68
CHAPTER 5 ....................................................................................................... 74
SOURCE LOCALIZATION ............................................................................................. 74
5.1 Introduction ............................................................................................................. 74
5.1. 1 Beamforming concept ......................................................................................... 75
5. 1.2 NearField Propagation Model ................................................... 77
5.2 NearField Propagation Model based Source Localization Algorithms .................. 80
5. 2. 1 Multiple Signal Classiﬁcation (MUSIC) ............................................................. 80
5.2.2 Conventional Focused Beamformer (CFB) ............................ 81
5.2.3 Linearly Constrained Minimum Variance Beamformer (LCMV) .................... 82
5.3 Algorithm Parameters .............................................................................................. 84
5.4 Algorithm Testing Procedure .................................................................................. 89
5.5 Results using the NearField Propagation Model based Techniques ...................... 92
5.5.1 Homogeneous Thorax ......................................................................................... 92
5.5.2 Mitral Valve Source .............................................................................................. 96
5.5.3 Right Lung Source ............................................................................................. 101
5.5.4 Discussion ........................................................................................................... 105
5.6 Transfer Function based Source Localization Algorithms .................................... 110
5.6.1 Transfer Function — Source Waveform Energy (SWE) ................................ 111
5.6.2 Transfer Function — MUSIC .............................................................................. 112
5.6.3 Transfer Function — CFB ................................................................................... 113
5.6.4 Transfer Function  LCMV ............................................................................... 113
5.7 Results using Transfer Function based Techniques .............................................. 1 15
5.7.1 Mitral Valve Source ............................................................................................ 115
5.7.2 Right Lung Source ............................................................................................. 118
5.7.3 Left Lung Source ................................................................................................ 121
5.7.4 Discussion ........................................................................................................... 125
5.8 Methods to improve NearField Propagation Model based Techniques ............... 127
5.8.1 Distributed Source Model ................................................................................. 127
5.8.2 Preliminary results using Distributed Source Model Techniques ............... 128
5.8.2.1 Homogeneous Thorax .................................................................................... 128
5.8.2.2 Mitral Valve Source ....................................................................................... 132
5.9 Conclusion ............................................................................................................. 137
CHAPTER 6 ..................................................................................................... 1 38
SUMMARY AND CONCLUSIONS .............................................................................. 138
6.1 Summary ................................................................................................................ 138
6.2 Conclusions and Future Work ............................................................................... 140
REFERENCES ................................................................................................. 142
vii
LIST OF TABLES
Table 31 Material properties (density, volume compressibility, shear elasticity, volume
viscosity, shear viscosity) and longitudinal wave velocity c, of the tissues that use the
Voigt model of viscoelasticity in the thoracic cavity. ....................................................... 22
Table 32 Material properties (density, volume compressibility, volume viscosity, shear
elasticity, shear viscosity) of air surrounding the thoracic cavity. .................................... 23
Table 33 Material properties (density, volume compressibility, shear elasticity, volume
viscosity, shear viscosity), longitudinal wave velocity CF and Young’s modulus E of the
spinal cord and invertebral disc in the thoracic cavity. ..................................................... 23
Table 34 Material properties (density, volume compressibility, shear elasticity,
resistance coefﬁcient to normal and shear stresses), longitudinal wave velocity c,, and
Young’s modulus E of the ribs and scapula bone that use the Maxwell model of
viscoelasticity in the thoracic cavity. ................................................................................. 24
Table 51 Localization error e in cm versus SNR in dB obtained in all three source
localization algorithms using I] = 0, for both pressure and velocity input signals in a
homogeneous thoracic medium. ........................................................................................ 95
Table 52 Localization error a in cm versus SNR in dB obtained in all three source
localization algorithms using I] = O, for both pressure and velocity input signals for a
mitral valve source ........................................................................................................... 100
Table 53 Localization error a in cm versus SNR in dB obtained in all three source
localization algorithms using 17 = 1, for both pressure and velocity input signals for a
mitral valve source ........................................................................................................... 100
Table 54 Localization error e in cm versus SNR in dB obtained in all three source
localization algorithms using 17 = O, for both pressure and velocity input signals for a right
lung source ....................................................................................................................... 104
Table 55 Localization error a in cm versus SNR in dB obtained in all three source
localization algorithms using I] = 1, for both pressure and velocity input signals for a right
lung source ....................................................................................................................... 104
Table 56 Localization error a in cm versus SNR in dB obtained in all three source
localization algorithms using r] = O, for only pressure input signals in a homogeneous
thoracic medium. ............................................................................................................. 108
viii
Table 57 Localization error e in cm versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a mitral valve
source. .............................................................................................................................. 1 18
Table 58 Localization error e in cm versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a right lung
source. .............................................................................................................................. 121
Table 59 Localization error a in cm versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a left lung source.
......................................................................................................................................... 124
Table 510 Localization error a in cm versus SNR in dB obtained in all distributed source
based algorithms using 71 = O, for both pressure and velocity input signals in a
homogeneous thoracic medium. ...................................................................................... 131
Table 511 Localization error a in cm versus SNR in dB obtained in all distributed source
based algorithms using a = O, for both pressure and velocity input signals for a mitral
valve source. .................................................................................................................... 135
Table 512 Localization error e in cm versus SNR in dB obtained in all distributed source
based algorithms using r] = 1, for both pressure and velocity input signals for a mitral
valve source. .................................................................................................................... 135
ix
LIST OF FIGURES
Images in this thesis/dissertation are presented in color.
Figure 31 Photographic crosssectional slice of the human thorax acquired from the
Visible Human Project database [40], [41] ....................................................................... 20
Figure 32 (a) Crosssectional slice of human thorax partitioned into distinct anatomical
regions, each tissue colorcoded with a unique color (b) Corresponding colorcode key to
identify the different tissues .............................................................................................. 21
Figure 33 Staggered grid layout for the 2D velocitystress FDTD scheme indicating the
relative spatial and time indices at which the variables are computed .............................. 34
Figure 34 Geometry of the problem for validating the FDTD code. Black circles
represent the array of receiver locations. Black cross represents the source location. ...... 41
Figure 35 j; components of the forcing functions (a) Gaussian ﬁrst derivative forcing
pulse (b) Sinusoidal forcing pulse ..................................................................................... 43
Figure 36 Comparison of vx and V; components obtained via FDTD to their
correspOnding analytical counterparts for the 1 kHz Gaussian ﬁrst derivative forcing
pulse (all colorbar units are in m/s) .................................................................................. 44
Figure 37 Comparison of vJr and v2 components obtained via FDTD to their
corresponding analytical counterparts for the 1 kHz sinusoidal forcing pulse (all colorbar
units are in m/s) ................................................................................................................. 45
Figure 38 Comparison of numerical FDTD with analytical results for a 1 kHz Gaussian
ﬁrstderivative point forcing fz source located at origin in a 2D homogeneous medium (p
= 1000 kg/m3, cp = 1500 m/sec, c3 = 500 m/sec) — Normalized vz component of velocity at
11 receiver locations located at z = 0.8 m and x = (1 m to 1 m, in steps of 0.2 m) .......... 46
Figure 39 Elastic wave propagation results as presented by Collino et al. [5 8] 
Snapshots of absolute velocity over the entire homogeneous domain caused due to an
explosive source located at the top left corner of the domain ........................................... 47
Figure 310 FDTD simulation results for the elastic wave propagation problem deﬁned in
Collino et al. [5 8] — Snapshots of absolute velocity over the entire domain caused due to
an explosive source located at the top leﬁ corner of the domain ...................................... 48
Figure 311 Source (white cross) and receiver (white circles and triangles) locations in
the thoracic geometry. (The axes are in units of m). ......................................................... SO
Figure 312 A 300Hz Gaussian ﬁrst derivative source signal j; used in the FDTD
simulations (a) Time waveform of the forcing signal (b) Corresponding Fourier spectrum.
........................................................................................................................................... 52
Figure 313 Magnitude of sourcetoreceiver transfer functions (TF) for a source at the
mitral valve location — TF plots at receiver #s 4, 8, 12, 16 and 20 in the 0100 Hz
frequency range ................................................................................................................. 53
Figure 314 Magnitude of sourcetoreceiver transfer functions (TF) for a source at the
mitral valve location — TF plots at receiver #3 4, 8, 12, 16 and 20 in the 100400 Hz
frequency range ................................................................................................................. 54
Figure 315 Magnitude of sourcetoreceiver transfer functions (TF) for a source at the
mitral valve location — TF image at all receivers in the 0100 Hz frequency range ......... 55
Figure 316 Magnitude of sourcetoreceiver transfer functions (TF) for a source at the
mitral valve location  TF image at all receivers in the 100400 Hz frequency range ..... 56
Figure 317 Snapshots (at 3, 7, 15, 25 and 40 ms) of velocity vJr and v2 distributions in the
entire thorax as a sound wave caused due to a 300 Hz Gaussian ﬁrstderivative forcing
source propagates ﬁ'om the mitral valve location. (The axes on each of the subplots are in
units of m and the colorbars have units of m/s) ............................................................... 59
Figure 318 Normal component of the velocity v,, at receiver #15 due to a 300 Hz
Gaussian ﬁrstderivative forcing source at the mitral valve location under three cases —
normal thorax, ribless thorax and shearless thorax ......................................................... 61
Figure 319 Normal component of the velocity v,. at receiver #12 due to a 300 Hz
Gaussian ﬁrstderivative forcing source at the mitral valve location for three different
thoraxsizes(L,XLz)—35cmX25cm,30cmX200mand40cmX3Ocm ................. 63
Figure 41 Anatomical simulation of a pneumothoracic air cavity on a 2D crosssectional
slice of the human thorax ................................................................................................... 66
Figure 42 Left column of images show the crosssectional slice of the thorax for varying
sizes of the pneumothoracic air pocket in the lung. The labels on the left side indicate the
size of the air pocket. The right column of plots represent the corresponding amplitude of
the source(o)tosensor(x) transfer functions. ................................................................... 69
Figure 43 Comparison of the amplitudes of the sourcetosensor transfer functions
corresponding to varying sizes of the pneumothoracic air pocket .................................... 71
Figure 44 Comparison of amplitudes of difference of the transfer ﬁmctions of varying
degrees of pneumothoracic severity and that of a normal thorax ...................................... 72
Figure 51 Concept of beamforming ................................................................................. 76
xi
Figure 52 Nearﬁeld propagation model .......................................................................... 77
Figure 53 The outer white box corresponds to the xz search space used in the
algorithms. The inner white box shows the size of the ﬁnest cell used in the search. ...... 85
Figure 54 Velocity 0 search values .................................................................................. 86
Figure 55 Sensor index and locations of the 19 sensor sites used in the source
localization algorithms ...................................................................................................... 87
Figure 56 An illustration of noisy waveforms at 10 dB, OdB and ~10dB compared to the
corresponding noiseless normalized pressure p and velocity v,. signals at four different
receiver sites  #s 9, 25, 41 and 57 ..................................................................................... 90
Figure 57 Location of the true source in the homogeneous thoracic medium is indicated
by the white cross. The white circles correspond to the sensor locations used in the source
localization codes. ............................................................................................................. 93
Figure 58 Objective function images of all three source localization algorithms using 27 =
0 obtained on both pressure and velocity input signals at an SNR of 5 dB in a
homogeneous thoracic medium. Both x and z axes of all images are in m. White cross
indicates the true source location while the white circle corresponds to the estimated
source location. .................................................................................................................. 94
Figure 59 Plot of localization error a in m versus SNR in dB obtained in all three source
localization algorithms using 17 = 0, for both pressure and velocity input signals in a
homogeneous thoracic medium. ........................................................................................ 95
Figure 510 Location of the true source at the mitral valve in the thoracic medium is
indicated by the white cross. The white circles correspond to the sensor locations used in
the source localization codes. ............................................................................................ 97
Figure 511 Objective function images of all three source localization algorithms using r]
= 1 obtained on both pressure and velocity input signals at an SNR of 10 dB for a mitral
valve source. Both x and z axes of all images are in m. White cross indicates the true
source location while the white circle corresponds to the estimated source location. ...... 98
Figure 512 Plot of localization error a in m versus SNR in dB obtained in all three source
localization algorithms using I] = 0 (top row) and r] = 1 (bottom row), for both pressure
and velocity input signals for a mitral valve source. ......................................................... 99
Figure 513 Location of the true source in the right lung of the thorax is indicated by the
white cross. The white circles correspond to the sensor locations used in the source
localization codes. ........................................................................................................... 101
xii
Figure 514 Objective function images of all three source localization algorithms using 17
= 1 obtained on both pressure and velocity input signals at an SNR of 10 dB for a right
lung source. Both x and z axes of all images are in m. White cross indicates the true
source location while the white circle corresponds to the estimated source location. 102
Figure 515 Plot of localization error e in m versus SNR in dB obtained in all three source
localization algorithms using r] = 0 (top row) and r] = 1 (bottom row), for both pressure
and velocity input signals for a right lung source..  .......................... 103
Figure 516 Location of the true source in the homogeneous thoracic medium is indicated
by the white cross. The white circles correspond to the sensor locations used in the source
localization codes. ........................................................................................................... 106
Figure 517 Objective function images of all three source localization algorithms using 7]
= 0 obtained on only the pressure input signals at an SNR of 5 dB in a homogeneous
thoracic medium. Both x and z axes of all images are in m. White cross indicates the true
source location while the white circle corresponds to the estimated source location. 107
Figure 518 All source and sensor locations considered for the transfer functions ........ 111
Figure 519 Objective ﬁrnction images of all four source localization algorithms obtained
on both pressure and velocity input signals at an SNR of 5 dB for a mitral valve source.
Both x and z axes of all images are in m. White cross indicates the true source location
while the white circle corresponds to the estimated source location ............................... 116
Figure 520 Plot of localization error a in m versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a mitral valve
source. .............................................................................................................................. 1 17
Figure 521 Objective function images of all four source localization algorithms obtained
on both pressure and velocity input signals at an SNR of 10 dB for a right lung source.
Both x and z axes of all images are in m. White cross indicates the true source location
while the white circle corresponds to the estimated source location ............................... 119
Figure 522 Plot of localization error e in m versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a right lung
source. .............................................................................................................................. 120
Figure 523 Location of the true source in the left lung of the thorax is indicated by the
white cross. The white circles correspond to the sensor locations used in the source
localization codes. ........................................................................................................... 122
Figure 524 Objective ﬁmction images of all four source localization algorithms obtained
on both pressure and velocity input signals at an SNR of 10 dB for a left lung source.
Both x and z axes of all images are in m. White cross indicates the true source location
while the white circle corresponds to the estimated source location ............................... 123
xiii
Figure 525 Plot of localization error a in m versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input sigials for a left lung source.
......................................................................................................................................... 124
Figure 526 Distributed source model ............................................................................. 128
Figure 527 Objective function images of all three distributed source based algorithms
using I] = 0 obtained on both pressure and velocity input signals at an SNR of 5 dB in a
homogeneous thoracic medium. Both x and z axes of all images are in m. White cross
indicates the true source location while the black diamond corresponds to the estimated
source location. ................................................................................................................ 130
Figure 528 Plot of localization error e in m versus SNR in dB obtained in all distributed
source based algorithms using 17 = 0, for both pressure and velocity input signals in a
homogeneous thoracic medium. ...................................................................................... 131
Figure 529 Objective function images of all distributed source based algorithms using r]
= 0 obtained on both pressure and velocity input sigrals at an SNR of 5 dB for a mitral
valve source. Both x and z axes of all images are in m. White cross indicates the true
source location while the black diamond corresponds to the estimated source location. 133
Figure 530 Plot of localization error a in m versus SNR in dB obtained in all distributed
source based algorithms using r] = 0 (top row) and r] = 1 (bottom row), for both pressure
and velocity input sigrals for a mitral valve source. ....................................................... 134
xiv
M‘Hl'ﬁl'ﬂqh
”S?”
.11)
0
mm
NOMENCLATURE
density
damping factor
velocity vector
stress tensor
force per unit volume vector
strain components
constitutive parameters
Lamé constants
stress
strain
elastic modulus
viscosity
viscosity in Maxwell’s model
viscosity in Voigt’s model
particle velocity in 1 direction
normal stress component in 1 direction
shear stress component in x—z plane
resistance coefﬁcient in the i:]' direction
forcing functions in the 1 direction
coefﬁcient of volume compressibility
coefﬁcient of shear elasticity
coefﬁcient of volume viscosity
coefﬁcient of shear viscosity
XV
31
“I
«E‘sTm“
558..
CID
fo
(x0, 2 0)
0'
B
angular frequency
wavenumber
stretching function in the I direction
scaling factor of the stretching ﬁmction
loss factor of the stretching function
stressrate tensor
strainrate tensor
stressrate tensor components
strainrate tensor components
thickness of PML
scaling order of the PML loss factor proﬁle
theoretical reﬂection coefficient of PML
longitudinal wave velocity
center ﬁ'equency of forcing function
xz location coordinates of source
standard deviation of Gaussian pulse
amplitude of forcing pulse
peak magritude at center frequency
timeshift of Gaussian pulse
unity disc with center (x0, 20) and radius a
gid step in x direction
gid step in z direction
gid step in time
stressrate cumulative sum components
strainrate cumulative sum components
xvi
maximum longitudinal wave velocity
shear wave velocity
slowest shear wave speed
number of spatial sample points to track cs, min
thorax dimension in x direction
thorax dimension in z direction
number of gid cells in x direction
number of gid cells in z direction
sampling frequency
normal component of velocity
sourcetoreceiver transfer function
mean pressure
.th . .
1 sensor srgnal mput to beamformer
sensor sigial input vector
number of sensors
. .th
beamformer weight of the 1 sensor
beamformer weight vector
output sigral of beamformer
beamformer output vector
velocity of wave propagation
number of acoustic sources
. th
distance of k source from reference sensor
. th .th
distance of k source from 1 sensor
th
k source waveform
source waveform frequency
xvii
. .th
norse at the 1 sensor
nearﬁeld propagation model steering vector of km source
nearﬁeld propagation model steering matrix
geometric spreading loss factor
exponential loss factor
sensor data covarianCe matrix
th . . .
k ergen vector of data covariance matrix
nearﬁeld propagation model steering vector corresponding to the
location (r, 0, CI) in a 3D polar coordinate search space
near—ﬁeld propagation model steering vector corresponding to the
location (x, z) in a 2D cartesian coordinate search space
objective function of the source localization algorithms
constraint steering vector matrix
constraint values vector
rectangular search space dimension in x direction
rectangular search space dimension in z direction
sigraltonoise ratio
localization error
th h .
k SOUI'CCtOlt sensor transfer ﬁrnctron
th .
k sourcetosensor transfer function vector
th .
z component of the k source ﬁmctron
. th . . h .
estimate of k source function usrng f sensor s1 gral and
corresponding transfer function
number of closely spaced omnidirectional point sources
xviii
CHAPTER 1
INTRODUCTION
1.1 General
Chest auscultation offers vast amounts of information that can assist in the diagrosis
of various diseases of the heart and lungs [1], [2]. Speciﬁcally, a multisensor acoustic
chest pad could be an extremely lowcost, powerful, lightweight and portable tool for
accurate quantitative diagnosis of thoracic pathologies. Existing technologies like
echocardiogaphy are fairly expensive and are typically not available at all places, for
instance, in outpatient home monitoring or battleﬁeld situations.
Recent years have seen many studies that attempt to exploit the spatiotemporal
distribution of thoracic sounds on the human torso for characterizing different
pathologies in patients with heart or lung diseases. Most of these rely on pattern
recogiition techniques or array—signal processing algorithms or a combination of both.
Although there has been some success in these efforts, not many of these studies have
addressed the effects of heterogeneity of the medium on the performance of the
algorithms [3], [4]. In addition, at lower audible range ﬁequencies, the effects of shear
waves in the sigrals measured cannot be igrored [5], [6]. The complexity of the problem
increases further because of the extreme nearﬁeld conditions that prevail. The algorithms
would thus seem to beneﬁt if they incorporate a more accurate physicsbased model, one
that governs the sound propagation in the human thoracic cavity.
In order to have a better understanding of the transmission of sounds in the human
thoracic cavity, this thesis attempts to numerically model the sound propagation using a
ﬁnitedifference time—domain (F DTD) approach laying emphasis on the true geometry
and material properties of the various tissues and organs in the human torso. The model
thus allows one to study the effects of the heterogeneities on the acoustic sensor sigrals
and can be used to validate all the existing sigral processing algorithms used in previous
studies. It also allows one to develop and validate newer signal processing and inverse
problem solutions to localize the pathological sites in the thorax in a noninvasive
manner. In addition, it can also be used as a tool to simulate different pathological
conditions of the heart and lungs and allow one to study the associated sound distribution
around the torso.
1.2 Dissertation layout
The layout of this dissertation is as follows:
Chapter 2 discusses some of the prior studies done in the ﬁelds of thorax acoustic
modeling and signal processing used for diagnostic purposes. Chapter 3 provides details
of the development of the viscoelastic FDTD model for sound propagation in the human
thorax proposed in this work and discusses some of the modeling results obtained using
the model. Chapter 4 employs the model developed to study the pneumothorax condition
of the lung and provides some results that might be of diagnostic utility. Chapter 5
discusses some of the traditional source localization algoritlnms that are validated using
the simulated signals obtained via the FDTD model and compares their performance to a
transferﬁnnction based algorithm developed using the forward model. Concluding
remarks are provided in Chapter 6 that also discusses the future steps to be taken to
address some of the unresolved issues.
CHAPTER 2
PRIOR WORK
2.1 Introduction
Studies related to understanding the origin, transmission, and characteristics of
auscultation sounds in addition to its diagnostic utility have been going on for more than
a few decades now. There have been various models proposed that attempt to simulate
the sound wave propagation in the thorax, and a thorough survey of these works allow us
to have a better understanding of the issues involved in pursuing this task. In addition to
the modeling work, there are many studies in literature that have focused their efforts to
developing signal processing algorithms either for source localization or for improved
diagnostic ability.
2.2 Acoustic Modeling Studies of the Thorax
The earliest studies [7]  [9] lay emphasis on estimating the transfer function of the
heartthorax acoustic system and understanding the frequency response and attenuation
characteristics of the heart sounds. These studies revealed that the heart sounds, in the
low frequency range of 50 Hz to 400 Hz, are damped by both the viscosity and the
geometrical spread in the medium, and attenuated by around 30 dB — 40 dB as they
reached the chest surface. Phonocardiogarn recordings indicated that almost 80% of the
power of sounds recorded on the chest wall arose due to a linear transmission of sounds
from the left ventricle. The remaining 20% are due to thoracic noise, noncoherent
contributions from extraneous signals and ﬁ'om system nonlinearities.
Verburg’s study [10] was one of the most signiﬁcant works of the 1980s that
attempted to model the heart sound propagation using viscoelastic wave equations,
although in a homogeneous media. A dipole model of the heart sound source was
proposed and the analytically simulated results of the model compared well with
experimental observations in many aspects. The thoracic media was considered to be a
linear viscoelastic body with complex material properties to account for the viscosity of
the tissues. It was observed that both shear and compressional waves propagated. The
study concluded that the lung tissue and the muscular connective tissue of the thoracic
wall are the two primary tissues to be considered in the media for wave propagation
modeling. These results formed the basis of all future thorax modeling efforts.
While efforts to understand the transmission of heart sounds continued, studies
focused on lung sounds were beginning to gather more interest in the scientiﬁc and
medical community. Wodicka et al. [11] — [15] developed an acoustic sound transmission
model from the respiratory tract to the chest wall via an equivalent acoustic circuit,
wherein the circuit elements represented different transmission characteristics of the
parenchyrna and the chest wall. The resonance peaks of the transmission spectra could be
estimated fairly well using these models. In addition, the models also estimated the
ﬁ’equency dependent propagation time (or phase delay) of sonic transmission from the
trachea to the chest wall. However, all of these studies assume a very simple
axisymmetric cylindrical geometry of the torso. Also, the tissue pathway to the chest is
treated simply as a mass load on the parenchyma.
Vovk et a1. [16]  [17] employed a similar axisymmetric cylindrical model to model
the acoustic properties of the chest in order to study the trarnsmission of breath sounds to
the chest. His work however, did incorporate the tissue pathway to the chest as a layered
structure of muscular rib cage, muscular fatty tissue and the skin. The study’s elegance
lay in the fact that for each layer, an appropriate set of differential equations of motion
was analytically solved to characterize the behavior of that layer. Also, the effect of
mechanical loading of the pickup transducers on the measured surface sounds was
careﬁrlly analyzed. It was observed that the presence of a transducer inevitably leads to
the conversion of longitudinal wave into surface waves. This can pose signiﬁcant
challenges when multisensor signal processing techniques are employed for analysis.
Vermarien too, in [18] had observed that the mecharnical loading of the sensors on to the
chest wall distorted the recordings of the sensor. Hence, the loading effect of the pickup
sensors should be taken into account in any experimental study.
In an effort to employ sound measurements to differentiate between normal and
diseased lungs, Leung and Sehati [19], estimated the average speed of propagation of
sound ﬁ'om the trachea to various sites on the chest wall over a range of frequencies
based on a crosscorrelation function of the input and the measured sound. The speed
estimates over a region of the tlnorax has been proposed to be used for diagnostic
purposes.
A one dimensional model for the propagation of pressure waves through the lung
using an equivalent mass spring chain model assuming a bi—periodic stack of tissue and
air layers has been studied in [20]. The work revealed that while the lung does show
dispersive behavior, i.e. signals of different ﬁequencies travel at signiﬁcantly different
speeds, for sufficiently slow rising pressure waves, the lung could well be approximated
using an equivalent homogeneous material in modeling studies.
The last few years have seen signiﬁcant developments in the ﬁeld of thorax modeling,
in particular, the study of sound propagation through lungs. Royston et al. [5] have
developed a 1D axisymmetric model similar to the work described in [16], accounting for
both attenuation and viscous damping losses in the medium to study the effect of
pneumothorax (PTX) on sound propagation ﬁ'om the bronchial airways to the skin
surface. In addition, a ﬁnite element model simulating a nonaxisymmetric 2D geometry
was used to study the same problem. The complex 2D geometry resulted in increased
nonradial scattering of sound than predicted by the 1D simple axisymmetric theory, thus
resulting in increased attenuation.
Royston et al. in [21] points out that at low audible frequencies, compression and
shear wave propagation from point sources can both be signiﬁcant although the shear
component becomes less signiﬁcant relative to the compression wave component beyond
a few hundred Hz. Royston et al. [22] extended their work in [5], [21] by modeling a lung
phantom and the true lung geometry (obtained from an Xray CT image) using a
boundary element formulation. Experimental results were compared with modeling
estimates for the lung phantom. It was observed that the shear waves were prominent
especially in the nonparenchymal soft tissue regions composed of fat, muscle and
visceral material. Also, the acoustic response ﬁeld created by an internal acoustic source
was signiﬁcantly different ﬁom that predicted by a ﬁ'eeﬁeld propagation model, thereby
emphasizing the importance of incorporating realistic geometries for modeling. Furtlner,
the study also investigated the possibility of using acoustic sound meaSurements on the
parenchymal surface to localize the lung sound source.
The ﬁnite difference time domain model developed by Narasimhan et al. [23] used
the actual tlnoracic geometry of the thorax obtained ﬁom the Visible Human Project to
model the sound propagation in the thorax. However, the model assumes a ﬂuid medium
throughout and ignores elastic wave contributions for the sake of simplicity. Key
observations of this work are:
0 Resonance effects of the thorax are signiﬁcantly lesser at lower frequencies (~
100 Hz) and are prominent at higher excitation frequencies (> 500 Hz). Also,
frequencies of the order of 100 Hz are most effectively transmitted through the
thorax.
o The size of the thorax has a profound effect on the transmission of breath sounds.
0 The spectrum of the human thoracic recordings is far richer relative to the
spectrum predicted by a homogeneous model simulation. This is probably due to
the various anatorrnical features in the thorax.
0 Chest wall sounds are affected by sound generated at the airways as well as the
manner in which they are propagated in the thorax. Hence, accurate models for
both sound generation and propagation are needed to simulate breath sounds.
Banks and Luke [24] employed the quasilinear viscoelastic equation proposed by
Fung [25] in a 2D ringlike geometry to study the propagation of shear waves in
biotissue. The model can be used in developing a noninvasive method to diagnose
coronary artery diseases caused by stenosis. The study justiﬁes the use of a viscoelastic
model as opposed to an elastic model for modeling shear wave propagation in biotissue.
All of the above discussed prior acoustic modeling studies, while, provide signiﬁcant
insight into the characteristics of sound propagation in the thorax, typically they either
assume a simple geometry or ignore the shear wave contributions [12], [16], [22], [23].
The work proposed in this dissertation attempts to overcome these limitations by
numerically modeling the viscoelastic sound propagation using a ﬁnitedifference time
domain (F DTD) approach laying emphasis on the true anatomical geometry and
experimentally obtained estimates of viscoelastic material properties of the various
tissues and organs in the human torso.
2.3 Signal Processing Methods for Acoustic Diagnostics
One of the earliest studies [26] on multisensor chest auscultation involved the
simultaneous recording of precordial vibrations on 22 sites on the chest and displaying
the peak amplitude spatial distribution as a function of time, which was used to analyze
the radiation and propagation patterns of the thoracic sounds. However, as mentioned in
[27], signal processing techniques that make use of the underlying physics, measurement
dynamics and noise models would prove to be more useful in practice. Emphasis needs to
be on modelbased signal processing strategies that can provide more meaningful and
useful results.
One such model based approach discussed in [3] uses the Multiple Signal
Classiﬁcation (MUSIC) algorithm for localization of acoustic sources in the heart by
employing a nearﬁeld propagation model. The main drawbacks of this approach are that
it assumes:
o a homogeneous medium of propagation
0 there is no scattering or reﬂection of sound in the thoracic medium
The authors extend their work in [28] employing the same model, but use a Choi
Williams distribution to compute a spatial timefrequency matrix that provides important
clues on the likely source locations in the heart.
Kompis et al. in [29], [30] and [31] employ a similar simpliﬁed model assuming
homogeneous wave propagation and a constant damping factor in their leastsquares
estimation approach.
Owsley in [4], [32] uses multiple auscultation recordings for localizing the artery
blockage sites. A nearﬁeld focused beamformer is used to perform imaging of the spatial
10
shear wave energy ﬁeld intensity. The beamformer model, as in [3] assumed a spherical
wavefront impingement along with homogeneous propagation characteristics.
To summarize, all the existing signal processing strategies for source localization in
the thorax [3], [4], [26] — [35] assume a constant propagation velocity in their models.
Although some of the algorithms do provide marginally decent results, signiﬁcant
improvement would be necessary before any of these techniques could be employed for
diagnostic purposes in a clinical environment. The viscoelastic forward model of sound
propagation in the thorax proposed in this dissertation, being based on wave propagation
physics speciﬁcally designed for human thorax, could thus potentially form the basis of
alternate signal processing algorithms for source localization and source waveform
estimation that might potentially yield better results as compared to traditional
algoritlnrns, eventually leading to clinical diagnostic solutions.
11
CHAPTER 3
VISCOELASTIC FDTD MODEL FOR HUMAN THORAX
3.1 Introduction
Auscultation sounds typically contain most of their energy at the lower audible range
of frequencies. At these frequencies, it is observed that compression and shear wave
propagation can both be signiﬁcant in tissue like materials. In general as the frequency
increases beyond a few hundreds of Hz, the shear wave component gadually becomes
less important. However, unlike the case of ultrasonic imaging, at the lower range of
frequencies, one cannot necessarily assume that compression waves dominate the
response over the shear and surface waves. In addition, the acoustic paths are far more
complex at sonic frequencies due to multiple reﬂections and standing wave patterns
within the thoracic cavity. Thus, in order to develop a reasonably accurate sound
propagation model in the human thorax, one is required to take into account botln the bulk
and shear modulus of elasticity and viscosity along with realistic estimates of the
geometries of the various tissues and anatomical structures in the thorax.
3.2 Governing Equations
As mentioned in the previous section, there is sufﬁcient evidence of the transmission
of shear waves in addition to longitudinal waves in the human torso at the lower audible
range ﬁ'equencies. Hence, the governing equations of sound transmission in the torso
requires to incorporate both these types of waves, and this is achieved through the
12
Elastodynamic equations, otherwise krnown as Elastic Wave equations. In addition, the
viscosities of the tissues are to be accounted for into these equations by suitably
modifying the coefficients or adding appropriate terms.
3.2.1 Elastodynamic Equations
These equations draw their roots ﬁom the Navier Stokes equation (force equation)
and the Hooke’s Law of elasticity [36], [37]. The Navier Stokes equation in its nonlinear
form is as follows:
(1)
Here p , I] , i3 , T and j7 denote the density, damping factor, velocity vector, stress
tensor, and the force per unit volume vector, respectively. (Note: For ﬂuid ﬂow, equation
(1) would include terms with shear and bulk viscosity. Also, tissues show viscoelastic
properties. One of the ways to incorporate damping characteristic of the viscous medium
into the above equations is through the damping factor 7) ).
In 2D, the equations can be written as:
6v, av, av, 6T,“ 5TH
—+ —+ — + =—+——+
p at vx 6x V2 52] 77vx 6x 62 fx
6v, avz av, 5TH 6TH
—+ —+ — + ——+
)0 at Vx V2 62 J 2 6x 5‘2 fz
(2)
The effect of the nonlinear term (I7.V)i3 was independently studied at the audible range
of frequencies as a 1D problem. The study revealed that its effect was negligible and the
13
additional level of complexity introduced due to the nonlinearity may well be ignored for
all practical purposes.
The Hooke’s law provides the relationship between the stresses T and strains S
acting on an element through the constitutive parameters. In its simplest form, for an
isotropic medium, the constitutive parameters that describe the medium completely are
functions of the Lamé constants (11, ,u).
(3)
Ignoring the nonlinear terms in equations (2) and combining the equations (2) and (3), we
can write the elastic wave equations for a 2D problem in terms of the velocity and stress
variablesas:
6v, 6TH 6TB
—+ + +
pa: nvx ax 62 fx
6v, 6sz 6T2,
+ _
pa: 2 ax 62 f2
(4)
arm av
—= 2+2 —£+2—z
,, ( ,) 0,,
6T2, 6v 6v
=2—x+ 2+2 —2
a, 0,, ( y)
at 62 6x
(5)
Equations in (4) are called the momentum equations, while equations in (5) are termed as
the constitutive equations.
14
3.2.2 Popular Viscoelastic Models
Although the equations in (4) and (5) could be utilized for modeling the viscoelastic
wave propagation, the parameters corresponding to the damping factor 7] aren’t available
for most tissues in the thorax. A standard linear viscoelastic, as opposed to the
elastodynamic model, is thus more appropriate for this problem as it incorporates the
viscous loss in the tissues in addition to their density and elasticity properties via
parameters that are more widely available in literature [38]. The Voigt and Maxwell’s
models are two of the most commonly employed rheological models of viscoelasticity of
tissues for numerical modeling purposes [25]. While the former is the more widely used,
the latter is also used for certain tissue structures like bones [39].
Both tlnese models use springs and dashpots to simulate the elastic and viscous
components of the stressstrain (CS8) relation. The spring component, associated with the
elasticity of the medium obeys the following relation:
0=E£
(6)
where E corresponds to the elastic modulus of the medium. The dashpot component,
associated with the viscosity of the medium 77 is governed by:
(7)
The Maxwell model consists of a spring and dashpot in series. The governing stress
stress equation for this model thus becomes:
15
60'
_1_ a
Etfi
5‘35
6:
(8)
However, the Voigt model consists of the spring and dashpot in parallel, and the
corresponding stressstrain equation thus becomes:
~66
0=E£+ —
”at
(9)
Since most viscoelastic studies on tissues employ either one of these two models to
incorporate the viscosity parameter, it would beneﬁt to combine these two models into
one framework that allows one to use parameters of either model. If we deﬁne ﬁ M and
ﬂy as the viscosity parameters associated with Maxwell and Voigt models respectively,
then the governing equations can be written as:
60' E as ,.. 62.9
—+:,——0' = E—+nV —
at ”M at azt
(10)
3.2.3 Viscoelastic Wave Propagation Equation
In order to incorporate the viscosity parameters associated with either of the Maxwell
and Voigt models into the governing equations, we employ the stressstrain relation as
described in equation (10) together with the momentum equations in (4) and obtain the
following equations that govern the wave propagation in a 2D viscoelastic medium:
16
+
p at ax 62 f"
av, 613,, 67,,
 —F
p at ax 62 f2
0TH 6v 6v
— ~+ 7‘ ==.t —41 + )1 +n2. —31
at 722 22 1 ax (1 #1) OZ
arm av av 23 av av
6t 7"”2 “[62 6x] ”zatiaz 6x)
(11)
(12)
(13)
(14)
(15)
where v, (I = x, z) are the particle velocities in the 1 direction, T]; are the normal stress
components and y” are the corresponding resistance coefﬁcients in the 1 direction, T m is
the shear stress and 1&2 is the corresponding resistance coefﬁcient in the xz plane, f; are
the forcing functions in the 1 direction, p is the density and 1:11 H (012, p=,u1+i 0112
are the complex Lamé constants of the medium. [I] and ,u] deﬁne the stiffness of the
17
isotropic medium under consideration arnd are termed as coefﬁcients of volume
compressibility and shear elasticity, respectively. ’12 and [12 deﬁne the viscosity of the
medium as pararneterized by the Voigt model and are termed as coefﬁcients of volume
viscosity and shear viscosity, respectively. a) is used to represent the angular frequency
of the wave. The resistance coefﬁcients of stress, originating ﬁ'om the Maxwell’s model
of viscosity, also control the viscous absorptive loss in the medium and are related to the
attenuation coefﬁcients of the tissues. Due to the heterogeneity of the medium, the
density p, Lamé constants (k and p) and the resistance coefﬁcients yij are all functions of
space governed by the viscoelastic properties of the various tissues in the thoracic cavity.
Equations (11) and (12) are called the momentum equations, while equations (13)  (15)
are termed as the constitutive equations.
18
3.3 Thorax Anatomy and Material Properties
In order to obtain an accurate structural representation of the various anatomical
features of the thorax, a photogaphic image of the thoracic cavity of a 39 year old male
human cadaver was acquired from the United States National Library of Medicine’s
Visible Human Project database [40], [41]. A trarnsverse crosssectional slice at the mid
thorax level passing through the heart, as shown in Figure 31, is chosen as the 2D image
for the model study. The various anatomical features, inclusive of the organ and tissue
boundaries are identiﬁed and the image is partitioned accordingly. Figure 32 shows this
partitioned colorcoded image where each color represents a distinct anatomical feature.
The material properties corresponding to these features are assigned to the respective
partitioned regions. The density, elasticity and viscosity parameters for each of these
tissues are obtained from existing literature [5], [12], [42] — [55] and are listed in Table 3
1 to Table 34. Some of these properties are observed to have a wide range of values
across different publications. In such cases, the most consistently reported value has been
chosen in the present work.
19
Right
Figure 31 Photographic cross—sectional slice of the human thorax acquired from the
Visible Human Project database [40], [41]
20
Air 
Ribs
Fashia — Subcutaneous Tissue
Senatus Anterior Muscle
Left Atrium
Erector Muscle. T Muscle
Pericardium
Atrium
Left Ventricle
Bone
lntervertebral Disc
Vein
Skin
Latissimus Dorsi Muscle
Cord
Pectoralis
Ventricle
Aorta
(b)
Figure 32 (a) Crosssectional slice of human thorax partitioned into distinct anatomical
regions, each tissue colorcoded with a unique color (b) Corresponding colorcode key to
identify the different tissues
21
Table 31 Material properties (density, volume compressibility, shear elasticity, volume
viscosity, shear viscosity) and longitudinal wave velocity cp of the tissues that use the
Voigt model of viscoelasticity in the thoracic cavity.
P a A c a
Material 3 1 ”I ”2 "
(kg/m) (Pa) (Pa) (Pas) (m/s)
Fashia 985 2.1146 5k b 4 ° 1465
Chest e e e
d 1050 2.6216 8.68k 9.73 1580
Muscles
Heart 1060 2.6336 124k f 15.9m f 1576
Blood 1060 2.6606 0 3,2m g 1584
Skin/ h 1
Dennis 1090 2.839G 2M 0,131 1615
Lungs 250 ’ 124.3k 41¢J 20.02 1‘ 23 ’
112 = 0 Pas [48]
ap and 6,, are obtained from Reference 46 unless explicitly mentioned
Reference 47
0Reference 48
dIncludes Serratus Anterior, Erector Spinae, Trapezius, Latissimus Dorsi and
Pectoralis Major Muscles
6Reference 49
{Reference 50
gReference 52
[Reference 53
Tefamce 12
J Reference 54. Shear wave velocity of 4 m/s is used to calculate the shear elasticity.
kReference 5
22
Table 32 Material properties (density, volume compressibility, volume viscosity, shear
elasticity, shear viscosity) of air surrounding the thoracic cavity.
Material P 3 11 12
(kg/m ) (Pa) (Pas)
Aira 1.24 147k 0.13
#1 = #2 = 0 Pa8 [42]
aReference 42
Table 33 Material properties (density, volume compressibility, shear elasticity, volume
viscosity, shear viscosity), longitudinal wave velocity cp and Young’s modulus E of the
spinal cord and invertebral disc in the thoracic cavity.
2. c E
Material p 3 I ”I "
(kg/m ) (Pa) (Pa) (In/S) (Pa)
Spinal Cord 1038 a 2.4676 466.7k 1542 ’ 1.4M b
Intervertebral c c
Disc 1000 16.4M 335.6k  1M
12 = 0 Pa—s [48]
pg = 4 Pas [48]
aReference 46
bReference 55
cReference 51 — Poisson ratio of 0.49 is also provided in this reference.
23
Table 34 Material properties (density, volume compressibility, shear elasticity,
resistance coefﬁcient to normal and shear stresses), longitudinal wave velocity cp and
Young’s modulus E of the ribs and scapula bone that use the Maxwell model of
viscoelasticity in the thoracic cavity.
P 2'1 [‘1 cp d E
(kg/m3) (Pa) (Pa) (m/s) (Pa)
Material
Ribs, Scapula 1990” 14.816 5.0636 3540a 13.96b
7n: yzz = 1.57k Us [45]
yxz =10.1k1/s[45]
21Reference 43
bReference 44
24
3.4 Boundary Conditions
The boundary conditions play a key role any numerical modeling scheme. In the
present application of viscoelastic wave transmission in the human thorax, boundary
conditions need to be imposed at the chestair interface. There are numerous ways to
incorporate appropriate boundary conditions at this interface. The most common
approach involves imposing a stressfree condition on the torso surface. However, for a
ﬁnitedifference approach using a rectangular gid, imposing a stressﬂee condition on a
complex torso boundary often becomes challenging. Instead, we propose to use a
perfectly matched layer (PML) of air surrounding a thin layer of air around the torso,
which can be incorporated into the ﬁnitedifference framework seamlessly.
3.4.1 Concept of PML
The basic idea of PML is to create an additional layer of a reﬂectionless absorbing
material around the geometry function that absorbs the waves in this region. This is
achieved through the use of complex stretching coordinates for the spatial variables. To
see how this works, let us consider a plane wave solution in a one dimensional x
coordinate space given by (“arena where a) and k represent the angular frequency and
wavenumber, respectively. This represents a nonattenuating plane wave. If the space
~ . w
variable x is replaced by a complex stretched version of it, say x = axx = [1+ri]x ,
a)
tlnen the plane wave solution modiﬁes to e—iwteikxe—kwxx/a). The third term
e_kwxx/ a; in this solution represents a loss factor that attenuates the wave with
25
increasing distance. Thus, by choosing a nonzero wx in the PML layer and setting the
remaining areas to have wx = 0 , we can simulate an absorbing material at the exterior.
It has been shown that an optimum choice of the stretching factors 8x can yield the
desired reﬂectionless property of the PML region [56].
Most PML formulations employing this coordinate stretching concept however,
require an artiﬁcial splitting of velocity and stress variables in the FDTD implementation
that decreases the computational efﬁciency of the codes. In order to avoid the
computational burden caused due to splitting, a nonsplit PML based on recursive
integation has been employed in the current work [57].
3.4.2 Recursive Integration PML
In the frequency domain, the 2D stretching functions 81(1 = x, 2) used in the PML can
be written as:
1W1
£1: a, +
(16)
where the real part a] is the scaling factor and the imaginary part w; corresponds to the
loss in the PML region. 6) represents the angular frequency. Correspondingly, the spatial
derivatives in the stretchedcoordinate PML space are given by:
2:. => —1—(: (I = x,z)
61 6‘1 31
(17)
26
The tilde here indicates that the spatial derivative is taken in the ﬁ'equency domain. The
recursive integation PML technique employs two auxiliary tensors, the stressrate tensor
S and the strainrate tensor E through which the nonsplit set of equations are obtained.
These tensors are deﬁned as follows:
.3 1 6T~ . .
U =—'a_fj (1,} =x,z)
Sj _]
(18)
Eij _1—§V—.l (i,j=x,z)
6'}: 6]
(19)
where the tilde corresponds to the frequency domain counterparts of the corresponding
variables. By introducing these tensor deﬁnitions and using equations (16) — (19) into the
frequency domain counterparts of equations (11) and (13), we get the following
equations:
iprx = Sxx + sz +fx
(20)
inxx +7):me =(/i.1+2/11)Exx +11Ezz
+ia)(/I.2 +2ﬂ2)§xx +i012§zz
(21)
Transforming these back into timedomain, and following the same procedure for all the
equations in (11) — (15), we obtain the following velocitystress equations:
av
p'a_tx=Sxx+sz+fx
(22)
27
6v
2 =Szx'l'Szz'ffz
6t
(23)
OT,“
“Et— “1' 7xxTxx=(/11 '1‘ 2.111 We: '1' AlEzz
6E 6E
+ 2 +2 —— 2 ——Z§
(2 #2) a 2 a:
(24)
OT
6:2 +7zszz = AlExx +011 4'21”] )Ezz
6E 6E
+2 —"1‘—+ 2 +2 J
2 at (2 412),),
(25)
6TH 6E 6E
—+ T Exz +E + ”+ 7*"
at 7:62 xz=#1( 2x) #2[ at at ]
(26)
The auxiliary tensor components here are obtained by substituting the deﬁnition of the
complex stretching ﬁrnction equation (16) into equations (18) — (19) and transforming
them into timedomain. The following equations result after some rearrangement:
wj LSU (z')dr=a; S,j (i,j=x,z)
(27)
wj [Eij (r)dr=:v—J Eij (i,j=x,z)
(28)
28
The recursive property of the integals become apparent following the discretization of
the above equations and is presented in section (3.6.1).
Equations (22) — (26) and (27) — (28) thus effectively incorporate the PML boundary
into the original viscoelastic equations. The values of a] and w] in the entire domain
decide the PML and nonPML regions. A value of a1 = 1 and w) = O is chosen in the
computational nonPML region, and a positive w; in the PML region to introduce the
loss. The reﬂectionless property of the PML is theoretically achieved only in a
continuous spatial coordinate system. While employed in a FDTD scheme, where the
domain is discretized, zero reﬂection cannot be practically achieved. In order to reduce
the numerical reﬂections in the PML region, the loss factor W] in this region is usually
chosen to have a tapered proﬁle as follows:
w, (u) = %10g[%)(%) (0 S u S d)
(29)
where u = 0 corresponds to the location of the interface of the computational domain and
the PML region, (1 represents the thickness of the PML, m is the scaling order, R is
the theoretical reﬂection coefﬁcient and CI, is the longitudinal wave velocity in the PML
region [56], [58], [59]. The amount of reﬂection can be controlled by choosing an
optimal R value. In the application presented in this paper, R and m were chosen to be
0.001 and 1 respectively [58].
29
3.5 Source Model and Initial Conditions
Sound generation in the thorax is an extremely complex phenomenon involving
various sources each of which is governed by different mechanisms. For instance, the
sounds originating at the heart alone are composed of the opening and closing beats of
the four valves, the vibrations of the ventricular and vascular walls and the ones caused
due to blood ﬂow haemodynannics. Similarly, modeling the respiratory breath sounds is
equally or more challenging as it involves complex air ﬂow dynamics in the trachea and
bronchial airways. In the current work, since the primary focus has been on sound
propagation rather than generation, the sound source model employed is a simple point
source placed at the mitral valve location, say at (x0, 20). A ﬁrst derivative of a Gaussian
pulse, centered at frequency f0 is used as the z component of the forcing function f; while
the z component of the forcing function ﬁr is set to zero.
P (tto)2 1
1 2
J2na
fz(x,z,t)=B(t—t0) e 20 5(XXO)§(ZZO) VtZO
=0, Vt <0
fx(x,z,t) = O, Vt
(30)
where the standard deviation 0' and amplitude B are deﬁned such that the peak magnitude
at the center ﬁequencyfo in the Fourier domain is kept constant. The forcing function f;
in the ﬁ'equency domain is given by:
30
—a'2a)2 . 7!
N 2 2 ﬂkao +3)
fz(x,z,a))=BO' we e 6(x—x0)6(zzo)
(31)
Thus, by deﬁning the standard deviation 0' and amplitude B as follows,
1 1
0' = — = ——
(00 279’ 0
(32)
B = A e05
0'
(33)
a peak magnitude of A is established at the center frequency 1?), regardless of the
frequency f0 The timeshift to is heuristically chosen to be 5 times the standard deviation
0' to provide a causal signal so that forcing function approaches to zero at t= 0.
Alternatively, instead of a point source, an explosive source may be modeled to
simulate a distributed source as follows:
76.1) = f(t)8'(r)
_ (tto)21
f(t)=B(t—to) 1 e 20’ 86—51086—2047120
27rd
=0,_Vt<0 _
(34)
31
where 13(0) is a urnity disc with center (x0, 20) and radius a.
The initial conditions for the problem deﬁne the velocity and stress variables at time t
= 0. Since the forcing function as deﬁned in equation (30) is zero for time t < 0, we
assume that the medium is in equilibrium at time t = 0, i.e. all the velocity and stress
variables are set to zero at t = 0 everywhere in the medium. Consequently, all the
components of the auxiliary tensors are also set to zero at t = 0.
32
3.6 FiniteDifference TimeDomain Implementation
There are a slew of ﬁnitedifference timedomain (F DTD) schemes available in
literature. However, the most commonly employed technique for elastodynamic wave
propagation problems is the center differencing scheme with a staggered gid that results
in a leapﬁog system of equations for the velocity and stress variables [60]. This provides
a second order accurate solution to the problem. In order to guarantee stability of the
FDTD codes, care must be taken in choosing the smallest timestep while handling the
discontinuities of material properties at the interface boundaries.
3.6.1 Discretization scheme
A staggered gid leapfrog centerdifferencing scheme [61] was used to discretize the
momentum and constitutive relation equations described in (22)  (26). For this 2D
velocitystress formulation of the viscoelastic equations, the velocity and stress variables
are used in a staggered gid as shown in Figure 33. In the discretized equations to follow,
the indices i, j and k correspond to the discretization of the spatial variables x, z and the
time variable t, respectively. Ax and AZ are the gid steps for the xaxis and Zaxis
respectively. At is the gid step irn time. Due to the leapfrog scheme employed, the
velocity variables and the stress variables are computed at times separated by half a time
313’.
ep2.
33
Z AH A F . vx
141» if? ‘ if 72? Txx’ Tzz
A T
1 A r XZ
1' 1+1
x———)
k k+1—
2
r A \f A \
J Txx' Tzz Sxx’ sz Exx' Ezz vx J
. 1 . 1
J +5 Szx’ Szz sz Vz Exz’ sz J +3
1 . l
l l+ 1 1+—
2 2
Figure 33 Staggered grid layout for the 2D velocitystress FDTD scheme indicating the
relative spatial and time indices at which the variables are computed
Using the center differencing scheme for each of the temporal and spatial derivatives
in equations (22) — (26), we get a sequence that constitutes the leapfrog system of
equations. Discretized forms of equations (22) and (24) are shown below:
34
. 1 . 1 . l , 1
Vx(l +§,j,k+§) =Vx[l +§,j,k—§)
+At Sn(i+%,j,k)+sz[i+—:,j,k)
1
+ ‘+—, ',k
fx(I 2 J )
_ 

b I
(35)
1313:0213")
. . 1 '
l (2.1+2pl)Exx(r,j,k+§)
1 7n .. I)
— —— +2E , ,k+—
[At+2)_ 1226.] 2 _
'31; [. k+l)_4e 111]
xx 2], 2 xx 2], 2
(ialﬂjzm +Exx(i,j,k—%]
3B [1’ 'k+l]—4E (i 'k—l)
22 3.], 2 22 a], 2
l 2At . . 3
(At +133] +Ezz[z,j,k—2)
_I
(36)
Note that the second terms on the left hand side of equations (24) — (26), i.e. the
Maxwell’s viscosity terms involving the stress variables are discretized by taking an
average of their values at t = km and t = (k+1)At as the stress terms are not evaluated at
(k+1/2)At. This averaging has the same secondorder accuracy as the other terms
involving the center differencing, and thus does not degrade the overall accuracy.
Similarly, time derivatives of the strainrate tensor components, i.e. the Voigt’s viscosity
35
terms in equations (24) — (26) are discretized using a threepoint stencil to provide
secondorder accuracy. 1
The integal of the auxiliary tensors in equations (27) — (28) are discretized using a
trapezoid rule and appear as follows for the Sn and En components.
(4+Wat’IS(..:,,,kl=I’”(‘+W>Tnl
2 Ax
wx . 1 .
—S +—, ,0
2 ”G 2’)
‘A‘ kl 1
+wJr Z Sxx(i+E,j,m)
m=l _
(37)
—v i+1 ‘k+1— 1
Wm .. 1 1 x 2’1’ 2
ax+ En r,1,k+— =—
2 Ax x14} .1.)
v ,j,k+—
_I. 2 2 ..
"xw '1
7Exx ”(1., j: '15]
Atk_1 1
, +
”if “(“m 2)
(33)
By imposing the initial conditions on the stressrate and strainrate components, and by
deﬁning the stressrate sums and strainrate sums as follows,
Qxx(i+:—,j,k1)= 19,211,,“ [1+],13m J
. 1 . . 1 .
=Q:(ll+§,j,k2J+Wx5xx[l+2,_],k1)
(39)
36
1 k“ 1
‘Pxx(i9j9m +5] = wx Z Exx(i:jam+§]
m=l
. . 1 . . 1
= xx('a.lam“2')+wxExx(la./ak_§]
(40)
we can rewrite equations (37) — (3 8) as
“+wa $304,130:T”('+1’J’k)—T’°‘(””k)
2 2 Ax
—Atﬂxx[i+%,j,kl)
. 1 . . l . . 1 .
Qxx l+§,_],k =0,“ l+§,j,k—1 +waxx r+—2,J,k
(41)
*v i+l 'k+— 
4 wat) 4 1) 1 x 2”’ 2
ax+ En 1,},k+— =—
2 2 Ax — (illelJ
_ x 2, , 2 ..
At‘I’xx(i,j,m+%]
‘I’ [i 'm+3—)‘P (i 'm+l)+wE [i 'k+—1)
xx 9.]: 2 ﬂ 9,], 2 x xx 3.], 2
(42)
Equations (41) — (42) provide a convenient way to recursively compute the integals
involved in the update equations of Sxx and Em. All the rernairning components of the
auxiliary tensors are discretized and updated in a similar fashion. The above equations
thus provide a timestepping system to track the propagation of viscoelastic waves
throughout the medium initiated by the forcing function.
37
3.6.2 Numerical Stability
For the above described explicit 2D FDTD scheme, one of the stability criteria is
provided by the CourantLevy condition that imposes a limit on the smallest timestep
allowed. According to the Courant condition,
cpAt\/ 1 2 + l 2 $1
(Ax) (A2)
(43)
must be satisﬁed everywhere in the medium. Thus, if Ax = A2 and cam is the
maximum longitudinal wave velocity in the medium, then the Courant numerical stability
condition reduces to
msi Ax
ﬁcp,max
(44)
Equation (44) is thus used to decide on the timestep for the FDTD scheme.
However, the Courant limit on the smallest timestep need not be sufﬁcient to
guarantee stability in a heterogeneous medium. Schroder [62] shows that the smallest
timestep required for a stable solution in elastodynamic wave propagation problems can
be more restrictive than the Courant condition as deﬁned in equation (44) depending on
the Lamé constants and densities of the materials at the interface regions in the medium.
To ensure that equation (44) is sufﬁcient for stability in any medium, the material
properties are averaged at the interface regions, especially at boundaries where the
properties differ signiﬁcantly. While the material densities are averaged linearly in x and
38
2 directions for the momentum equations, the inverse of the shear elasticity p1 of each of
the four adjoining cells is averaged for use in the constitutive equations.
39
3.7 Numerical Results
The leapfrog system of FDTD equations developed in section (3.6.1) is implemented
in Matlab. The codes are validated on a simpliﬁed geometry before using them for
simulating the sound propagation in the thorax from a source placed at the location of the
mitral valve in the heart. Source to receiver transfer ﬁrnctions are computed and analyzed
to provide some perspective on how the results may be used for diagnostic purpose.
3.7.1 Validation of FDTD codes
In order to validate the FDTD simulation codes, the numerical results obtained are
compared to analytical expressions derived for a simpliﬁed case. We consider a
homogeneous, inﬁnitely extending, isotropic, and elastic medium with a density of 1000
kg/m3, and longitudinal and shear velocities of 1500 m/s and 500 m/s with a point force
source acting at the origin. Two types of forcing functions are considered and the
comparison results with both of these shall be presented here. In addition to comparison
with analytical solutions, the FDTD results, particularly in terms of the effectiveness of
its PML layers, are compared with those presented by Collino et al. in [58].
3.7.1.1 Comparison with analytical solution
The two fz forcing functions used to compare the FDTD results with analytical
solutions are:
0 First derivative of a Gaussian pulse
0 Sinusoidal pulse
40
In both cases, the j; component of the forcing function is set to zero. The velocity and
stress components are calculated at an array of 21 receiver locations at z = 0.8 m, x = [1,
~09, O.8,. . .,0.1, O, 0.1, ..., 0.9, l] m in the x—z coordinate system as shown in Figure 3
4.
Geometry: Homogeneous medium 3
cp=1500 rn/s, cs=500 m/s, p=1000 kg/m
1.5 r . . . .
1  J
ooooooooooooooooooooo
0.5  
E
N
0 r x 
0.5  
x (m)
Figure 34 Geometry of the problem for validating the F DTD code. Black circles
represent the array of receiver locations. Black cross represents the source location.
The 2D FDTD program discretized the domain into 120 X 120 cells with a PML
width of 18 cells on each side of the domain. A cell size of 16.66 mm X 16.66 mm and a
41
timestep of At = 1 us is used in the simulation. The displacement solutions for the elastic
wave equations for an impulsive point force as derived by Eaton et al. are used to obtain
analytical results for the applied forcing functions, which are then compared with the
FDTD simulation results [63]. Figure 35(a) shows a ﬁrst derivative of a Gaussian pulse
with a center frequency of 1 kHz that is used as the f; component of the forcing function
for the ﬁrst case while Figure 35(b) shows the 1 kHz sinusoidal pulse f; forcing function
used in the second case. The velocity components vx and vz obtained using the FDTD
simulations for the two forcing functions are compared as images in Figure 36 and
Figure 37 [64]. Figure 36 compares the results obtained for the Gaussian ﬁrst derivative
forcing pulse while Figure 37 compares those for the sinusoidal forcing pulse.
Figure 38 compares the normalized vz component obtained using the analytical
expression and the FDTD simulation for the 1 kHz Gaussian ﬁrst derivative forcing pulse
at 11 of the receiver locations, corresponding to x = [1, 0.8, 0.6, 0.4, 0.2, 0, 0.2, 0.4,
0.6, 0.8, 1] m, by overlapping the 1D signal waveforms obtained in both cases in the
same plot. We observe excellent agreement between the numerical and analytical results.
The largest normalized mean squared error (NMSE) across all sensors was 2.5% and the
mean NMSE was 1.55%. Further, we also notice that the PML is effective in absorbing
the outgoing waves as no reﬂections are observable in Figure 38 [58].
42
Force (N/ma)
0 1 2 3 4 5
Time (ms)
(8)
I;
2
A 1
g
E,
0
2
O
u.
I
—l
rt:
0 2 4 6 8 10 12
Time (ms)
0))
Figure 35]; components of the forcing functions (a) Gaussian ﬁrst derivative forcing
pulse (b) Sinusoidal forcing pulse
43
Receiver
er#
Receiv
vxFDID vzFDTD
. 0.6
#15 4‘ = 0.4
. .5 0.2
‘ o
.02
' o.4
1 2 3 4 5 1 2 3 4 5
Time (me) Time (ms)
vx  Analytical vz  Analytical
. 0.6
0.4
f” 0.2
O
0.2
' o.4
1 2 3 4 5 1 2 3 4 5
Time (me) Time (ms)
Figure 36 Comparison of Vx and v2 components obtained via FDTD to their
corresponding analytical counterparts for the 1 kHz Gaussian ﬁrst derivative forcing
pulse (all colorbar units are in m/s)
vxFDTD vzFDTD
20
2
E15
3 o
310
0
ﬂ: 5 _2
1 2 3 4 5 1 2 3 4 5
Time (ms) Time (ms)
vx Analytical 122Analytical
20
2
#15
5
.2
81
¢
:1:
20
4:15
S
0 '310
d)
n:
.2 5
1 2 3 4 5
2 3 4 5
Time (ms) Time (ms)
Figure 3 7 Comparison of Vx and v2 components obtained via F DTD to their
corresponding analytical counterparts for the 1 kHz sinusoidal forcing pulse (all colorbar
units are in m/s)
45
Normalized vz
— Numerical
  Analytical
11 —
10
31: 9 
8 
0) 6—\/\—’/\t
4 M
3 M
2 2/\_\/\,
1 —A— \/\r
0 1 2 3 4 5
Time (ms)
Figure 38 Comparison of numerical FDTD with analytical results for a 1 kHz Gaussian
firstderivative point forcingﬁz source located at origin in a 2D homogeneous medium (,0
= 1000 kg/m3, cp = 1500 m/sec, cs = 500 m/sec) — Normalized vz component of velocity
at 11 receiver locations located at z = 0.8 m and x = (l m to l m, in steps of 0.2 m)
3.7.1.2 Comparison with results published in literature
Francis Collino et al. presented some results on the application of PML to linear
elastodynamic problems in [58]. The propagation of elastic waves in a homogeneous
elastic medium due to an explosive source at the left top comer is shown as a series of
time snapshots of the absolute velocity over the entire domain in Figure 39 and Figure 3
46
10; the former one being the result published in [58] and the latter one, obtained here
using the previously described FDTD code.
t = 2.76ms t = 5.53ms t = 8.351118
t: 11.1m3 t = 13.8ms t : 16.71718
Figure 39 Elastic wave propagation results as presented by Collino et al. [58] —
Snapshots of absolute velocity over the entire homogeneous domain caused due to an
explosive source located at the top leﬁ comer of the domain
47
t = 16.7ma
Figure 310 FDTD simulation results for the elastic wave propagation problem deﬁned in
Collino et al. [58] — Snapshots of absolute velocity over the entire domain caused due to
an explosive source located at the top left comer of the domain
3.7.2 Simulation results on human thorax
Having validated the FDTD codes for stability and accuracy, they are then used to
model the sound wave propagation in the thoracic cavity by discretizing a transverse
crosssectional slice image of the thorax acquired from the Visible Human Project as
described in section (3.3). The spatial grid size was determined based on the number of
spatial sample points (K S) chosen to track the shortest wavelength in the medium. Thus,
48
if cs, min is the slowest shear wave to be traced ande is the bandwidth of the forcing
function, then the grid size is obtained by
Ax=Az=1—Cs’min
KS fmax
(45)
The rectangular domain of Lx X L2 = 35 cm X 25 cm enclosing the thorax was thus
discretized into NJr X N2 cells using the grid size as determined in equation (45). PML
layers, 5 cm thick, surrounded this rectangular domain on each side. The timestep was
chosen to be 0.9 times the Courant limit as deﬁned in equation (44). The forcing function
fz was a ﬁrst derivative of a Gaussian pulse with a center ﬁ'equencyﬁ) placed at the mitral
valve location as described in section (3.5). The 2D F DTD codes were used to calculate
the velocity and stress components at an array of locations (sensor locations) along the
torso boundary and were stored at a desired sampling frequency F S Figure 311 shows
the locations of the source and receivers in relation to the thoracic geometry. The location
of the source is indicated by a white cross and the receivers on the front chest surface are
represented by white circles. All the remaining receiver locations are represented by
white triangles. The receivers are indexed in a clockwise manner, the index numbers of
few of which are indicated in the ﬁgure.
49
Figure 311 Source (white cross) and receiver (white circles and triangles) locations in
the thoracic geometry. (The axes are in units of m).
The simulations used a center frequency f0 of 300 Hz that yielded a 40 dB bandwidth
of around 1 kHz. The forcing ﬁmctionfz and its frequency spectrum are shown in Figures
312(a) and 312(b) respectively. The sourcetoreceiver transfer ﬁinctions (TF) were
computed for each of the receiver signals. The normal component of the velocity v” was
used for the receiver signals. Figure 313 and Figure 314 show the magnitude of the
transfer functions of receiver #5 4, 8, 12. 16 and 20 in the frequency range of 0100 Hz
50
and 100400 Hz, respectively. Figure 315 and Figure 316 show the same for all the
receivers around the torso as an image. The following parameter values were used in the
FDTD codes to obtain these results [65].
Grid size: Ax = A2 = 1mm
Dimensions : N x = 350, N z = 250, PMLwidth =150 layers
Time step: At = 0.15,usec
Source ﬁmction : f0 = 300Hz, A =106
Sampling frequency : F S = 8000Hz
(46)
51
x 10"
5  4
«6"
E
\
a 0
RN
5 4 _
0 10 20 30 40
Time (ms)
o 200 400 600 800 1000 1200
Frequency (Hz)
(b)
Figure 312 A 300Hz Gaussian ﬁrst derivative source signal]; used in the FDTD
simulations (a) Time waveform of the forcing signal (b) Corresponding Fourier spectrum.
52
Le—R4
.5
(a)
N
VeocityForce TF(00)
A
I
20 4'0 60 80 00
Frequency (Hz)
Figure 313 Magnitude of sourcetoreceiver transfer functions (TF) for a source at the
mitral valve location — TF plots at receiver #5 4, 8, 12, 16 and 20 in the 0100 Hz
frequency range
53
—l
1
.0
m
.9
m
p
a.
VeocityForce TF(<0)
0.2 ~
100 150 200 250 300 350 400
Frequency (Hz)
Figure 314 Magnitude of sourcetoreceiver transfer functions (TF) for a source at the
mitral valve location — TF plots at receiver #5 4, 8, 12, 16 and 20 in the 100400 Hz
frequency range
54
VeocityForce TF(co)
20
11?
E 40
>.
O
C
d)
g. 60
9
U.
80
100
20 40 60
Receiver#
Figure 315 Magnitude of sourcetoreceiver transfer functions (TF) for a source at the
mitral valve location — TF image at all receivers in the 0100 Hz frequency range
55
VeocityForce TF(o))
200
250
Frequency (Hz)
300
350
10 20 30 40 50 60 70
Receiver#
Figure 316 Magnitude of sourcetoreceiver transfer functions (TF) for a source at the
mitral valve location — TF image at all receivers in the 100400 Hz frequency range
From Figure 313 and Figure 314, it is clear that the heart~tochest acoustic transfer
functions at different receiver locations along the front surface of the thorax differ
signiﬁcantly, although, all of them have most of their energy content lying in the O80Hz
frequency range. In Figure 314 we notice that all the receivers except #R8 have peaks
56
between 125 Hz and 170 Hz, and a distinct valley at around 180 Hz. Similar observation
can be made ﬁom the image in Figure 316 across all the 20 receivers on the ﬁont
surface. Also, receivers #8 to #15 have high levels of acoustic transmission in the 180
270 Hz range as compared to all other receiver locations. In Figure 3—15 we observe that
the 040 Hz signals transmit well to the posterior left corner of the thorax corresponding
to receiver location #32 to #40. These observations although cannot be conclusively used
for characterizing the transmission paths of an S1 sound that originates from the mitral
valve source, they do provide some limited insight into the auscultation sites around the
torso that might be better suited for study of certain types of abnormal heart sounds.
The source model described in section (3.5) could be modiﬁed to simulate lung sound
propagation by modeling explosive sources as suggested in equation (34) at the locations
where the secondary bronchi meet the plane of crosssectional slice in the lungs. Since
the frequency range of lung sounds extend well beyond 1 kHz, the center ﬁ'equency
would then have to be set to a higher value of around 1 kHz to incorporate a wider
bandwidth into the forcing function [1]. Consequently, the time step and spatial grid size
would get reduced according to equations (44) — (45).
3.7.3 Wave Propagation Characteristics
In order to gain a better insight into the propagation of the viscoelastic waves
originating from the source, we track the velocity components in the entire domain as a
function of time. This allows us to observe the complex paths that the sound waves
traverse as they reach the surface of the torso. The effect of the elasticity and viscosity of
the various tissues in the thoracic cavity on the sound propagation can thus be studied
using this model in an effective manner. Figure 317 shows the vx and vZ component
57
distributions in the thorax at ﬁve time snapshots as the Gaussian ﬁrstderivative pulse
source waveform propagates ﬁom the mitral valve location towards the periphery of the
thorax [65]. At the 3 ms snapshot of the vz image, we see that the initial wavefront of the
sound waves have already reached the front surface of the torso right above the heart
valve. This speedy conduction can be attributed to the absence of the lung lobes in these
paths. The waves passing through the lungs are slowed down and severely attenuated as
shown in the vx and vz images of all the remaining time snapshots. The images at 15 ms,
25 ms and 40 ms also show effects of multiple reﬂections at the heartlung interface. A
detailed quantitative and qualitative analysis of these snapshots involving tracing the
waveforms in each tissue individually may be performed, if necessary.
58
0.1 1 4
E o " o 2
M
0.1 ,1 0
0 0.1
0.1 ' 0 5 0.1 , . 1
w r ‘ ' ‘
E o \’ U0 0 (Q /’ '05
(a . 0
0.1 0.5 0.1 ' 0.5
0 0.1 0 1 0 0 1
w 0.1 __ 0.4 01 (g ‘ 3.3
E , , 3;. .
2 0 0 0 0
0.1 ~04 0.1 .0.4
0 0.1 0 1 0 0 1
In 0 V: 0 0 0
N
0.1 _0_2 0.1 0.2
0 0.1 0 1 0 0 1
_ 1 0.1 
2 0.1 t . 0.05 0.1 0.05
5 o . 0 0 0
V 0.1 _ 0.05 0.1 0.05
0 0.1 0.1 0 0.1
Figure 317 Snapshots (at 3, 7, 15, 25 and 40 ms) of velocity vx and vz distributions in
the entire thorax as a sound wave caused due to a 300 Hz Gaussian ﬁrstderivative
forcing source propagates from the mitral valve location. (The axes on each of the
subplots are in units of m and the colorbars have units of m/s)
59
3.8 Effects of Anatomical Variations
The viscoelastic FDTD model of the thorax can be used to study the effects of rib
structures, shear elasticity of tissues, and size of the thoracic cavity on the signals
reaching the surface of the thorax by making subtle modiﬁcations to the parameters. The
study, whilst aids understanding the effects of such anatomical variations, also justiﬁes
the need to incorporate shear elasticity and rib structures into the thoracic geometry.
3.8.1 Effects of ribs and shear elasticity of tissues
It is clear from Table 31 to Table 34 that the ribs and bones have distinctly different
material properties as compared to the other tissues. One thus expects their presence to
cause a signiﬁcant impact on the sound waves reaching the chest surface. A simulation
study to observe the effects of ribs on the signals is conducted by assigning the material
properties of the muscles to the rib and bone structures. The waveforms thus obtained are
compared to the signals obtained on a normal thorax. Similarly, a simulation was carried
out to study the effect of shear elasticity of the tissues on the sound signals. This was
performed primarily to justify the inclusion of shear elasticity into the model. The
coefﬁcient of shear elasticity ,u] was set to zero in the entire domain to simulate the
shearless case. The normal component of the velocity V" obtained at receiver #15 for a
normal thorax, ribless thorax and shearless thorax are compared in Figure 318 [65].
The simulations used the same parameters as listed in section (3.7 .2) in equation (46).
Results indicate that the presence of ribs clearly attenuates the signals signiﬁcantly,
although the overall shape of the waveform is not affected as much. However, the
60
absence of shear elasticity in the thorax results in a completely different waveform shape,
thereby illustrating the importance of inclusion of shear elasticity in the model.
Since the shearwave contribution in certain tissues like the lungparenchyma are
insigniﬁcant compared to that of the compressional waves, more efﬁcient numerical
schemes could ideally be employed in such regions [5]. However, this would require a
more diligent analysis to tackle the internal boundary conditions. The FDTD numerical
scheme used in the current work simpliﬁes the task and provides a generic framework,
which may be applied to modeling various anatomical and pathological conditions too.
I 1
[ — Norma]
0'05 e— No Ribs
4 No Shear
0.05
Time (ms)
Figure 318 Normal component of the velocity v” at receiver #15 due to a 300 Hz
Gaussian ﬁrstderivative forcing source at the mitral valve location under three cases 
normal thorax, ribless thorax and shearless thorax
61
3.8.2 Effect of size of thorax
All the simulations in the preceding sections have assumed a ﬁxed rectangular
domain of Lx X L2 = 35 cm X 25 cm enclosing the thorax. As one would expect, the size
of the thorax might also have an effect of the receiver signals. In order to study this, the
original thorax image was linearly resized to ﬁt into 30 cm x 20 cm and 40 cm x 30 cm
rectangular domains and F DTD simulations were carried out with the same spatial and
time discretization. The normal component of the velocity v,, obtained at receiver #12 for
the three different thorax sizes are compared in Figure 319 [65]. Results indicate that
while the amplitude of the signal corresponding to a larger sized thorax is smaller and
timedelayed compared to that of the smaller sized thorax, the shape of the waveform is
retained atleast upto the ﬁrst waveﬁont, corresponding to the dominant longitudinal wave
peak and the shear wave peak. The secondary reﬂections caused by, presumably the rib
structures, are not as prominent in the 40 cm X 30 cm thorax as in the 30 cm X 20 cm
thorax.
62
0.04 I .
—35 cm x 25 Em
e—30 cm X 20 cm
A~40 cm X 30 cm
0.03
I
0.02
I
0.01
I
(m/S)
I
>= 0.01
0.02
I
0.03
I
0.04
I
l
0.05 ' 1 1
Time (ms)
Figure 319 Normal component of the velocity v" at receiver #12 due to a 300 Hz
Gaussian ﬁrstderivative forcing source at the mitral valve location for three different
thorax sizes (Lx X L2) — 35 cm X 25 cm, 30 cm X 20 cm and 40 cm X 30 cm
63
CHAPTER 4
STUDY OF SOUNDS UNDER PNEUMOTHORAX CONDITIONS
4.1 Introduction
Pneumothorax is a pathological condition of the lungs in which air collects in the
pleural space, i.e. the space that surrounds the lungs, leading to a lungcollapse
eventually. Since it is a potentially life threatening disease, and is easily treatable when
diagnosed early, there is signiﬁcant interest in developing a noninvasive, lowcost, quick
and advanced diagnostic scheme. Currently, the goldstandard technique for diagnosis of
pneumothorax is the chest Xray or a chest CT although misdiagnosis rates of around
30% have been reported even with these methods [66], [67].
Auscultation sounds offer a viable alternative for the diagnosis of various thoracic
pathologies including pneumothorax [68]. Prior experimental studies have shown the
effectiveness of breath sounds in differentiating between normal and pneumothorax
conditions. By digitizing these sounds using acoustic sensors mounted around the torso,
and combined ‘with an intelligent multisensor processing scheme, one could devise a
lowcost, noninvasive and portable solution for accurate diagnosis that may be even used
for outpatient home monitoring purposes. As mentioned in section (1.1) an accurate
sound propagation forward model simulating pneumothorax conditions can aid in the
development of such multisensor processing schemes.
The model would facilitate the study of effects of pneumothorax on the acoustic
sensor signals and can be used to develop and validate various multisensor processing
64
algorithms for diagnostic purposes. It also allows one to develop newer inverse problem
solutions to localize and quantify the size of the pneumothoracic air cavity in the thorax.
Section (4.2) provides simulations details pertaining to modeling the pneumothorax
conditions using the FDTD codes. Preliminary simulation results demonstrating the effect
of varying sizes of air cavity on the sourcetosensor transfer functions are presented in
section (4.3). Concluding remarks discussing the future steps to be taken to address some
of the unresolved issues are also discussed in section (4.3).
65
4.2 Simulation of Pneumothorax
As described earlier in section (4.1), a pneumothorax is basically an air pocket
trapped in the pleural space between the lung and the chest wall. The condition is thus
modeled by introducing a region at the exterior edge of the lung that has material
properties of air, as shown in Figure 41. In order to study the effect of the size of the air
pocket on the chest signals, sizes varying from 1.7 cm to 5 cm was introduced, size being
deﬁned as the maximum extent of the air cavity.
Pneumothoracic air pocket
Right side of the thorax : z axis (m)
0.15 0.1 0.05 0 0.05 0.1 0.15
Back side of the thorax : x axis (m)
Figure 41 Anatomical simulation of a pneumothoracic air cavity on a 2D crosssectional
slice of the human thorax
66
To model the sound wave propagation in the thoracic cavity, we ﬁrst discretize the
transverse crosssectional slice image of the thorax. The rectangular domain of Lx X L2 =
35 cm X 25 cm enclosing the thorax is discretized into Nx X NZ cells using the grid size
as determined in equation (45). PML layers, 15 cm thick, surrounded this rectangular
domain on each side. The time step was chosen to be 0.8 times the Courant limit as
deﬁned in equation (44).
The 2D FDTD codes were used to calculate the velocity and stress components at an
array of sensor locations placed all around the torso. Subsequently, the mean pressure, p
was computed using the normal stress components as follows:
p=—(Tn +Tzz)/2
(47)
The FDTD simulations assumed the Maxwell’s model of viscoelasticity for the entire
domain and used the following parameter values [69]:
Grid size :Ax = A2 = 1mm
Dimensions : N x = 350, N z = 250,
PML width = 150 layers
Time step :At = 0.15psec
Source function :fo = 300Hz, A =106
Sampling frequency :F S = 8000 Hz
7 values 37x2: = yzz = 7x2 = 300 Vx,z
(48)
67
4.3 Results and Conclusion
Using the simulated pressure waveforms at the sensor locations around the torso and
the source forcing function, the sourcetosensor transfer functions (TF). were computed
for each pneumothorax simulation. For ease of analysis, we consider only one of the key
sensor locations on the front surface of the chest, indicated by a blue cross in the left
column of thorax anatomy images of Figure 42. The source location is indicated by the
white circle in the images. These images represent the various pneumothorax conditions
simulated with increasing degree of severity going from top to bottom, the labels on the
left indicating the size of the air cavity. The right column plots in Figure 42 represent the
amplitudes of the sourcetosensor pressure TFs for the corresponding pneumothorax
condition. The image and TF plot in the ﬁrst row corresponds to a normal thorax. The
maximum amplitude of the TF obtained for the normal thorax is used as the scaling factor
for all the plots in this ﬁgure and subsequent ﬁgures.
68
Thorax Anatomy Normalized Pressure
_ 0.1 “‘12.. ". 1.W
E . , 3 ,
o 0.5
2
O 200 400
1
5
r~_ 0.5 
0 200 400
1
E
N 0.5 
N
0 200 400
1
g M
,_ (15 
co"
0 200 400
1
g MM
l0 2 V's, ..., , r
01 0 01 200 400
Both x & z axes in m Frequency (Hz)
Figure 42 Leﬁ column of images show the crosssectional slice of the thorax for varying
sizes of the pneumothoracic air pocket in the lung. The labels on the leﬁ side indicate the
size of the air pocket. The right column of plots represent the corresponding amplitude of
the source(o)tosensor(x) transfer functions.
69
In order to ﬁrrther study and compare these TFs, Figure 43 plots these on one
common axis. We notice that at frequencies around 50100 Hz there isn’t any signiﬁcant
difference in the amplitudes of the TFs. However, the frequency range of 350475 Hz
shows a distinct trend of decreasing amplitudes with increasing degree of severity of
pneumothorax. Figure 44 considers the amplitudes of the difference transfer functions
(DTF), computed as:
DT F = T F (Pneumothorax) — T F (Normal)
(49)
The plot shows that the deviation of the TF from the normal case is the largest for the
pneumothorax with the largest air cavity (5 cm, in this case), and shows a clear trend
across all the cavity sizes. More importantly, this deviation is enhanced in the 350475
Hz ﬁ'equency range, thereby providing some additional insight into possible diagnostic
utility of chest signals measured in this frequency range.
Similar results may be obtained from other sensor sites too. As one would expect,
sites located signiﬁcantly distant from the air cavity do not show such prominent
deviation. Thus, by combirning information from various sensor sites, one could devise a
multisensor noninvasive auscultatory system for diagnosis of pneumothorax conditions.
70
Normalized Pressure
+Normal
1 » +Pneumo 1.70m
+Pneumo 2.7cm
+Pneumo 3.70m
0.8 +Pneumo 5.00m
0.6
0.4
1" \ A
0.2 x" ‘
100 200 300 400 500
Frequency (Hz)
Figure 43 Comparison of the amplitudes of the sourcetosensor transfer functions
corresponding to varying sizes of the pneumothoracic air pocket
71
1 I l I I I
e Pneumo 1.70m
B Pneumo 2.70m
—e— Pneumo 3.70m
+ Pneumo 5.00m
.0
co
Normalized Pressure
.0 .0
h 0)
9
N
g
’1
0 H
100 200 300 400 500
Frequency (Hz)
.. 0 l‘
0
Figure 44 Comparison of amplitudes of difference of the transfer functions of varying
degees of pneumothoracic severity and that of a normal thorax
The 2D numerical FDTD viscoelastic model developed in chapter (3) has thus been
used to study the lung sound propagation in the human thoracic cavity under
pneumothorax conditions. The model provides a convenient testbed to investigate the
possibility of using an array of noninvasive acoustic sensors for diagnosing
pneumothorax and estimating the extent of severity of the disease. Future efforts would
be directed towards developing a 3D model of the torso that may provide a more accurate
72
representation of the thoracic sounds. Studies pertaining to effects of the location of the
cavity in the thorax would be considered. Furthermore, the models could be used to
design and evaluate various inversion schemes that can be used to estimate the location
and size of the air cavity.
73
CHAPTER 5
SOURCE LOCALIZATION
5.1 Introduction
Physicians and doctors listen to auscultation sounds and are able to diagnose various
heart and lung pathologies using single stethoscope measurements made at distinct sites
on the torso. This diagnostic ability can be signiﬁcantly enhanced vvitln a multisensor
system that makes simultaneous recordings at a bunch of auscultation sites around the
torso. In addition to allowing the physician to listen to or even view the sound signatures
of the simultaneously recorded sounds, a multisensor system can be used to combine the
sensor recordings in an optimal manner to localize the sound, i.e. identify the exact
locations of the intrathoracic acoustic sources. Given the locations of the sources, it may
then be used to amplify the sounds originating from only a particular source in the thorax
while suppressing the interfering sounds from other sources, thereby aiding physicians in
diagnosis. Source localization for thorax can thus assist in developing improved non
invasive diagnostic solutions.
The viscoelastic numerical model of the thorax discussed in the previous chapters can
be used as a convenient test bed for validating many of the traditional source localization
algorithms as well as development of newer and better schemes. The sections to follow
discuss the some of the algorithms used traditionally and assess their performance to this
application. Further, an alternate algorithm based on the viscoelastic model is proposed
74
and assessed. Finally, a modiﬁed version of the traditional algorithm is proposed to
alleviate some of the problems associated with it.
5.1.1 Beamfonning concept
The algorithms that combine the information from distinct, spatially separated,
simultaneous recording sensors to localize the sources are popularly termed as
bearnforrners [70]. The beamformer usually weights and combines the sensor signals in
an optimal fashion so as to extract the desired signal of a particular frequency originating
ﬁ'om a particular point in space while suppressing all contributions of signals from other
points in space. One of the challenges associated with beamfonning to a desired source
thus is, noise and interference from undesirable sources as illustrated in Figure 51. If
xi(n) denote the N sensor signals and Wi denote the corresponding weights, 1' = 1,2,. . .N,
then the beamformer output signal y(n) is given by:
N
y(")= éwz‘xiln)
(50)
The above is typically expressed in vector form as y = WHX, where w = [w 1, W2:
wN] is the weighting vector that is optimally determined by the beamforrning algorithms
to meet the desired response of the beamformer.
75
ARRAY OF SENSORS
L
L
DESIRED
SOURCE
Combine LOCATION
sensor “ESTIMATE
signals OF DESIRED
I optimally SOURCE
NOISE &
INTERFERING I
SOURCES
Figure 51 Concept of beamforrning
The basic element of all beamforming algorithms is a wave propagation model that
relates the source signal to the received sensor signals. Hence, the performance of the
beamformer not only depends on the strength of noise and interfering sources, but also on
relevance of the assumed wave propagation model. Many of the traditionally employed
beamforrning algorithms were originally developed for radar signal processing and far
ﬁeld applications where the wave propagation model essentially contains a phase delay.
However, the. current application of sound source localization in the thorax is
predominantly a nearﬁeld problem, which is far more complex. In recent times, there
have been attempts to modify the farﬁeld models to make them suitable to near ﬁeld
propagation problems [3], [4]. However, the extension from farﬁeld to nearﬁeld as
assumed by these algorithms hold good, if at all, only in a homogeneous medium of
transmission. Thus, in practice, the algorithms governed by the nearﬁeld propagation
76
model may have to be tweaked by controlling a few parameters that optimize their
performance in complex media of the thorax.
5.1.2 NearField Propagation Model
The nearfield propagation model assumed by the algorithms is based upon a
spherical wavefront transmission through a homogeneous medium of velocity C. Assume
there are N sensors receiving signals from K acoustic sources. Let r0k denote the distance
(In . th
of the k source from the reference sensor and rik denote the distance of the k source
.tl . .
from the l 1 sensor as shown below in F ngure 52.
sensor 1' reference
sensor
0 9 e
\O
O ‘. O
s. .
x. \
r,\e \rUk
1k \ .
‘. \
\ O
. \
\..
sourcek
Figure 52 Nearfield propagation model
If Sk(t) represents the km source waveform of frequency f, then according to the
model, the 1le sensor signal xi( t) can be expressed as in equations (51) and (52). The
77
noise contribution ni(t) at each sensor is usually modeled as a zeromean i.i.d
(independent identically distributed) Gaussian noise with variance 02.
K
xiii): ZSk(t)aik(rOkarikafat)+ni(t)a 1951"
(51)
where
,6
r r . r — r
aik (rOk ’rik ’f’t) = (i) exp[_ ﬂ—{LJCXP( 1271M)
rik r 0k 0
(52)
The vector ak = [a 1k, 02b a Nk]T is popularly termed as the steering vector, while the
aik’s describe the nearﬁeld propagation model [3], [71]. In equation (52), we notice that
the model contains three explicit terms that are described below:
(i) The ﬁrst term is related to the geometric spreading loss due to the spherical
wave expanding from an omnidirectional point source. The factor ﬂ, called the
geometric spreading loss factor takes a value of l for a homogeneous medium.
(ii) The second term is associated with the energy loss per unit volume in a lossy
homogeneous medium, and the factor n is correspondingly called the
exponential loss factor of the medium.
(iii) The last term corresponds to the phase correction associated with the time
delay for the signal to propagate ﬁ'om the source to the sensor. The average
velocity of propagation 6 plays an important role in this term.
78
Clearly, in the above model, the three factors (,3, 17, 0) do not possess a well deﬁned
meaning while being used in a complex heterogeneous thoracic medium. Hence these
factors are chosen from a range of values that optimize the performance of the
algorithms.
The following section shall describe three different source localization algorithms that
were implemented to validate their utility for localizing an acoustic source in the human
thorax. Results obtained using these algorithms have been discussed in the subsequent
section.
79
5.2 NearField Propagation Model based Source Localization
Algorithms
The three algorithms employing the nearﬁeld propagation model described in the
previous subsection considered are:
i) MUltiple SIgnal Classiﬁcation method (MUSIC)
ii) Conventional Focused Beamformer (CF B)
iii) Linearly Constrained Minimum Variance beamformer (LCMV)
5.2.1 Multiple Signal Classification (MUSIC)
The MUSIC algorithm [3], [72] is a subspacebased estimation technique that uses
the eigen structure of the sensor datacovariance matrix to determine the source locations.
Using equation (51), the data covariance matrix R = E{XXH} can be expressed as
follows:
R = AEth )1” + 021 = ARSAH + 021
(53)
where A = [a 1, (12, . . ., aK] is the steering matrix and I is the identity matrix obtained
assuming a spatially white Gaussian noise distribution. It is clear from equation (53) that
if there exists K incoherent sources, then the singular value decomposition of R will
result in K signal eigenvectors (with an eigenvalue ok2+02, k = l, 2,. . ., K where each
okz corresponds to the signal power of the kth source signal) and (N — K) noise eigen
80
vectors (each with an eigen value of oz). Thus, if vk is a noise eigenvector (i.e. k > K),
tlnen the following hold good:
Rvk = ozvk
H 2 .
=> (ARSA + o I)vk = o’zvk ﬁom equation (53)
H
=> ARSA V], = 0
=> AHvk = 0 (since Rs and A are botln of rank K)
This implies that the (N — K) noise eigenvectors are ortlnogonal to the K source
direction vectors, and lie in the nullspace of A. This forms the basis for most of the
eigenvector based algorithms. The MUSIC algorithm uses this concept and estimates the
source locations by identifying the peaks of a function deﬁned as follows:
1
N
vaa(r,t9,¢4
k=K +1
P(r, 6, ¢) =
2
(54)
where a( r, 0, D) is the nearﬁeld propagation model based steering vector
corresponding to the location (r, 6, [:1 ) in the search space. If a(r, 6, 1:] ) corresponds to
the steering vector from a source location, then the denorrninator of the function
approaches zero, thus peaking the function at that location.
5.2.2 Conventional Focused Beamformer (CFB)
The conventional focused beamformer, also called as Bartlett beamformer,
determines the weights deﬁned in equation (50) by maximizing the output power [72]. In
81
other words, the optimal weight vector W is obtained by maximizing the quantity
E {yyH}=WHRW. In order to obtain a nonuivial solution the weight vector is restricted
to have a norm of 1. Substituting the sensor signal model as deﬁned in equation (51) into
the output power expression, we can obtain the optimal weight vector W solution as:
a(r, 6, ¢)
Ja”(r,6.¢)a(r.6.¢)
w:
(55)
The expression for the output power thus becomes:
p(r,0,¢)  a” (r, 6’ ¢>Rao<”0000000“x
0 10 20 30 40
Velocity Index
Figure 54 Velocity c search values
The sensor locations play a key role in most source localization algorithms.
Optimization of sensor locations, however, is an extremely challenging problem and
beyond the scope of this work. Typically, most applications employ an evenly spaced
distribution of sensors to their problems. A similar approach has been chosen here.
Simulated signals ﬁ'om a set of 19 receiver sites, evenly spaced along the periphery of the
torso have been employed for source localization except for one case in section (5.5.1)
86
that uses only 5 sensors on the front surface of the torso. Figure 55 shows the receiver
index numbers and corresponding locations on the thorax surface.
Sensor Index and Locations
2 axis (m)
0.15 —0.1 0.05 0 0.05 0.1 0.15
xaxis (m)
Figure 5—5 Sensor index and locations of the 19 sensor sites used in the source
localization algorithms
In addition to the location of the sensor, the type of sensor measurement also plays an
important role in many applications. Pressure, velocity, acceleration, and displacement
are the most commonly measured quantities for sound studies. In the current work, both
the mean pressure p and normal component of the velocity v" are used for source
localization. Results obtained using both of these are reported.
87
Finally, the number of time snapshots and the timewindow of the sensor signals used
also dictate the performance of the algorithms. A timewindow of 16.5 ms sampled at 8
kHz was used in the current work. Other shorter windows were also evaluated; but none
of them yielded signiﬁcantly better results.
w"? .—
88
5.4 Algorithm Testing Procedure
The following test protocol was adhered to assess the performance of each of the
source localization algorithms discussed in this chapter.
0)
(ii)
(iii)
(W)
The sensor signals, botln p and v" were obtained from the viscoelastic FDTD
simulations that used a 300 Hz Gaussian ﬁrst derivative pulse as the f},
sourcing function.
White Gaussian noise of varying SignaltoNoise Ratio (SNR) levels, ranging
from 15 dB to 10 dB, was added to artiﬁcially introduce noise into the sensor
signals. Figure 56 compares the noisy waveforms of p and v" to the noiseless
signals obtained due to a nnitral valve source for 3 different SNR levels at
receiver #3 9, 25, 41 and 57.
The noisy signals are then ﬁltered using a 7th order lowpass Chebyshev ﬁlter
with the stopband ripple 20 dB down and with a cutoff ﬁ'equency of 1200
Hz to eliminate the high frequency contributions of the noise.
The noisy signals are then passed through the source localization modules.
Each of the source localization algorithms computes its objective function
over the xz search space for every 77 and velocity 6 as described in section
(5.3).
89
Normalized Pressure Normalized Velocity
57.  W
41   , 
25
Nolseless
9
57
41
10 dB
0dB
10 dB
0 10 0 15
Time (ms) Time (ms)
Figure 56 An illustration of noisy waveforms at 10 dB, 0dB and 10dB compared to the
corresponding noiseless normalized pressure p and velocity v" signals at four different
receiver sites  #s 9, 25, 41 and 57
90
(v) The objective function image is normalized and then smootlned using a 4 cm
X 4 cm averaging ﬁlter to eliminate any noisy peaks in the function.
(vi) The xz location of the maximum intensity in the smoothened objective
function for every I] and velocity 0 is stored. The median xz location across
the range of velocities is used as the ﬁnal estimate of the source location for
each I] and SNR level.
(vii) The localization error 8, i.e. the distance between the true source location and
the estimate is tabulated and plotted against the SNR values.
(viii) The above procedure from step (iii) to (vii) is performed for both pressure p
and velocity v" signals and for every source localization algorithm.
In the following sections, the results obtained using the MUSIC algorithm on pressure
p waveforms is titled as MUSICP while that obtained on velocity v" waveforms is titled
as MUSICV. Similar convention is adhered to for the remaining two localization
algorithms (CFB and LCMV).
91
5.5 Results using the NearField Propagation Model based
Techniques
This section presents the results of the MUSIC, CFB and LCMV algorithms
employing the nearﬁeld propagation model on the simulated signals for two different
source locations in the thorax. In addition, results on a homogeneous thorax are also
presented for comparison.
5.5.1 Homogeneous Thorax
A homogeneous thorax with a longitudinal velocity cp of 1605 m/s and density p of
1055 kg/m3 is considered. The medium is considered shearless and viscousless. In other
words, shear wave contributions are absent and the medium is lossless. Figure 57 shows
the locations of the true source and sensors used in the source localization codes.
92
Homogeneous Thorax
= 1605 m/s, c5 = 0 m/s,p—=1055 kg/m3
ooooo
O.1
A 0.05
.5,
.‘L’ 0
X
(6
N 0.05
0.1
0.2 0.1 0 0.1 0.2
x axis (m)
Figure 57 Location of the true source in the homogeneous thoracic medium is indicated
by the white cross. The white circles correspond to the sensor locations used in the source
localization codes.
Figure 58 shows the objective function images obtained for each of the source
localization algorithms using a loss factor 11 of 0, for both pressure and velocity input
signals at an SNR of 5 dB. The images in the ﬁgure also indicate the true source location
by a white cross and the estimated location by a white circle. Further, the title of each
image provides the localization error 8 of the corresponding algorithm.
93
MUSICP I a = 0.92195 cm MUSICV I e = 5.5902 cm
0.05 0.05
0 0
0.05 0.05
0.1 0.1
0.1 0 0.1 0.1 0 0.1
CFBP 1 e = 1.7117 cm
0.05
0
0.05
0.1 .
0.1 0 0.1 0.1 0 0.1
LCMVP I e = 0.92195 cm LCMVV Z a = 5.5036 cm
9
0.05 x
'0.1 0 0.1 01 0 0.1
Figure 58 Objective ﬁrnction images of all three source localization algorithms using 17 =
0 obtained on botln pressure and velocity input signals at an SNR of 5 dB in a
homogeneous thoracic medium. Both x and z axes of all images are in m. White cross
indicates the true source location while the white circle corresponds to the estimated
source location.
The localization error 8 of each of the algorithms using 77 = 0 is plotted against SNR
in Figure 59 for both the pressure and velocity input signals.
94
P I n = 0 V I r] = 0
e MUSIC 1» MUSIC
0.25 o CFB ~ 0.25. o CFB 
0.2, A LCMV 0.2 A LCMV
g 0.15 1* g 015
to 0 to
0.1’ 9 0.1 .
0.054 o  0.05tastes.
A A Q 9 9 O
0' . . . 1 0— . . .
10 0 10 10 0 1O
SNR SNR
Figure 59 Plot of localization error 8 in m versus SNR in dB obtained in all three source
localization algorithms using I] = 0, for both pressure and velocity input signals in a
homogeneous thoracic medium.
Table 51 presents the same results as illustrated in Figure 59 in a tabulated form.
Table 51 Localization error 8 in cm versus SNR in dB obtained in all three source
localization algorithms using I] = 0, for both pressure and velocity input signals in a
homogeneous thoracic medium.
SNR (dB)
15
10
0
5
10
I]
0
0
0
0
0
s MUSICP (cm)
1.97
1.97
0.92
1.71
15.16
8.19
s CFBP (cm)
2.69
2.69
1.71
4.94
12.59
7.83
s LCMVP (cm)
2.69
2.69
0.92
0.92
1.22
5.59
s MUSICV (cm)
5.59
5.59
5.59
5.50
8.81
5.59
a CFBV (cm)
5.50
5.50
5.50
5.53
8.81
3.45
a LCMVV (cm)
5.50
5.50
5.50
5.59
5.50
5.59
95
The above results indicate that the nearﬁeld model assumed by the algorithms ﬁt the
pressure signals better than the velocity signals. From Figure 58, it appears that although
the direction of the beamshape is fairly accurate, the range estimate obtained for the
velocity signals is erred. For the pressure signals, the source location estimates are
signiﬁcantly accurate, MUSIC and LCMV algorithms performing the best. All the
algorithms show a deterioration of performance with decrease in SNR level as expected.
5.5.2 Mitral Valve Source
A point source located at the nnitral valve in the thorax is considered. Figure 510
shows the locations of the true source and sensors used in the source localization codes.
Figure 511 shows the objective function images obtained for each of the source
localization algorithms using a loss factor I] of 1, for both pressure and velocity input
signals at an SNR of 10 dB. The images in the ﬁgure also indicate the true source
location by a white cross and the estimated location by a white circle. Further, the title of
each image provides the localization error 8 of the corresponding algorithm.
96
Mitral Valve Source
z axis (m)
0.15 0.1 0.5 0 0.05 0.1 0.15
xaxis (m)
Figure 5— 10 Location of the true source at the mitral valve in the thoracic medium is
indicated by the white cross. The white circles correspond to the sensor locations used in
the source localization codes.
97
MUSICP : a = 14.1156 cm MUSICV I e = 11.8004 cm
—0.05
0
0.05
0.1 .
0.1 O 0.1 0.1 0 0.1
CFBP : a = 4.1049 cm CFBV : e = 5.883 cm
O
0.1  .,
0.1 0 0.1 0.1 0 0.1
LCMVP I e = 7.7389 cm LCMVV I e = 11.0368 cm
30.1 0 0.1
Figure 511 Objective function images of all three source localization algorithms using 77
= 1 obtained on both pressure and velocity input signals at an SNR of 10 dB for a mitral
valve source. Both x and z axes of all images are in m. White cross indicates the true
source location while the white circle corresponds to the estimated source location.
The localization error 8 of each of the algorithms using 77 = 0 and I] = l is plotted
against SNR in Figure 512 for both the pressure and velocity input signals.
98
P I n = 0 V: n = 0
1r MUSIC 4 MUSIC
0.25 * O CFB 0.25 ' Q CFB
A
0.2 LCMV 1 0.2 A LCMV 1
$5015": n A. .s t 4' £045} A .s 5 4. 4‘
00 A to 4. A
0.1' o A  0.1
O O O O
0.05 O O O O O . 005 O n
0' . . J 0' . 0 I .
—10 0 10 10 0 10
SNR SNR
P I n =1 V I n =1
1r MUSIC 4. MUSIC
0.25 t O CFB 1 0.25  O CFB
A A
0.2 LCMV _ 0.2. LCMV .
w e e w 1' e 3 e
0.1 O A 0.1 A ‘
005 A 005 o I ° ° °
. ‘ O O . ' O ‘
O o 0 O
0* . . . .01 . . .
10 0 10 10 0 10
SNR SNR
Figure 512 Plot of localization error 8 in m versus SNR in dB obtained in all three
source localization algorithms using I] = 0 (top row) and 77 = 1 (bottom row), for both
pressure and velocity input signals for a mitral valve source.
Tables 52 and 53 present the same results as illustrated in Figure 512 in a tabulated
form.
99
‘_.:IiI.L
Table 52 Localization error 8 in cm versus SNR in dB obtained in all three source
localization algorithms using I] = 0, for both pressure and velocity input signals for a
mitral valve source.
SNR (dB)
15
10
5
10
'l
0
0
0
0
8 MUSICP (cm)
14.33
14.33
14.12
14.12
14.33
14.33
a CFBP (cm)
5.17
4.60
4.05
4.60
5.37
9.54
a LCMVP (cm)
10.82
9.84
14.12
14.12
14.12
14.12
a MUSICV (cm)
14.12
14.12
14.12
14.12
11.80
14.12
e CFBV (cm)
6.86
6.86
6.86
5.53
0.61
6.86
e LCMVV (cm)
14.33
14.12
12.06
14.12
14.12
14.12
Table 53 Localization error 8 in cm versus SNR in dB obtained in all three source
localization algorithms using I] = 1, for both pressure and velocity input signals for a
mitral valve source.
SNR (dB)
15
10
5
5
10
t]
1
1
l
l
1
e MUSICP (cm)
14.33
14.12
11.80
11.80
14.33
14.33
a CFBP (cm)
4.60
4.10
3.47
3.72
4.60
9.54
a LCMVP (cm)
9.27
7.74
14.12
14.12
14.12
14.12
r: MUSICV (cm)
11.80
11.80
11.80
12.57
7.00
14.12
e CFBV (cm)
5.88
5.88
6.58
4.74
2.51
6.58
e LCMVV (cm)
14.33
11.04
9.53
14.12
12.81
14.12
The above results indicate that the nearﬁeld model assumed by the algorithms result
in signiﬁcantly erroneous estimates of the source locations for both the pressure and
velocity signals. From Figure 511, it appears that neither the direction of the beamshape
100
nor the range estimates are accurate. Among the three algorithms, CFB provides the best
estimates, although not as accurate as desired. Also, the results appear to be marginally
better for r] = l as compared to 17 = 0.
5.5.3 Right Lung Source
A point source located in the right lung of the thorax is considered. Figure 513 shows
the locations of the true source and sensors used in the source localization codes.
Right Lung Source
“—5.. ,
5.:Tﬁ5'5, A" 
z axis (m)
0
I:
0.15 0.1 0.05 0 0.05 0.1 0.15
xaxis (m)
Figure 513 Location of the true source in the right lung of the thorax is indicated by the
white cross. The white circles correspond to the sensor locations used in the source
localization codes.
101
MUSICP 2 e = 6.5765 cm MUSICV I e = 7.1701 cm
0.05
0
0.05
. 0.1
0.1 0 0.1 0.1 0 0.1
CFBP : e = 8.6377 cm CFBV : e = 11.3864 cm
.'.‘ E
0.05 005
0 0
0.05 0.05
0.1 I 0.1
.0.1 o 0.1 0.1 0 0.1
Figure 514 Objective function images of all three source localization algorithms using 17
= 1 obtained on both pressure and velocity input signals at an SNR of 10 dB for a right
lung source. Both x and z axes of all images are in m. White cross indicates the true
source location while the white circle corresponds to the estimated source location.
Figure 514 shows the objective function images obtained for each of the source
localization algorithms using a loss factor 71 of 1, for both pressure and velocity input
signals at an SNR of 10 dB. The images in the ﬁgure also indicate the true source
102
location by a white cross and the estimated location by a white circle. Further, the title of
each image provides the localization error 8 of the corresponding algorithm.
P2n=0 V2n=0
0.3 . . 0.3 . .
t MUSIC 1* MUSIC
0.25 o CFB 0.25 o CFB
0.2, A LCMV , 0.2 A LCMV ,
€0.15 E015
w o 00 o 0 o o o
0'123 922 01445335”
A A A
0.05» 0.05
0»; . . og . .
10 0 10 10 0 10
SNR SNR
P1n=1 V2n=1
4 MUSIC 1* MUSIC
025' o CFB 025’ o CFB
0.2, A LCMV , 0.2, A LCMV ,
@015 £0.15
to to O O O O O
0.0513 A A Z A 0.05 8 A
0 . . . 0 . .
10 0 1o 10 0 10
SNR SNR
Figure 515 Plot of localization error 8 in m versus SNR in dB obtained in all three
source localization algorithms using I] = 0 (top row) and I] = 1 (bottom row), for both
pressure and velocity input signals for a right lung source.
103
The localization error 8 of each of the algorithms using 7] = 0 and I] = l is plotted
against SNR in Figure 515 for botln the pressure and velocity input signals. Tables 54
and 55 present the same results as illustrated in Figure 515 in a tabulated form.
Table 54 Localization error 8 in cm versus SNR in dB obtained in all three source
localization algorithms using I] = 0, for botln pressure and velocity input signals for a
right lung source.
SNR (dB) 15 10 5 0 5 10
,, 0 o 0 o o 0
S MUSICP (cm)
8.46
8.46
8.46
8.46
8.92
7.78
e CFBP (cm)
10.15
10.15
9.95
9.95
10.74
9.44
e LCMVP (cm)
6.61
6.03
6.58
8.46
8.30
7.71
c MUSICV (cm)
8.46
8.46
8.46
8.46
8.92
8.46
e CFBV (cm)
12.18
12.18
12.18
8.05
12.68
12.30
a LCMVV (cm)
8.30
7.80
7.16
8.46
8.92
8.46
Table 55 Localization error 8 in cm versus SNR in dB obtained in all three source
localization algorithms using 7] = l, for botln pressure and velocity input signals for a
right lung source.
SNR (dB)
15
10
5
10
r]
1
1
l
l
e MUSICP (cm)
6.58
6.58
6.58
8.46
8.92
6.61
e CFBP (cm)
8.39
8.64
9.17
9.17
10.02
8.73
s LCMVP (cm)
5.40
4.84
5.14
8.46
6.23
5.87
8 MUSICV (cm)
7.17
7.17
7.17
6.58
8.30
8.46
a CFBV (cm)
11.39
11.39
11.39
5.45
12.68
12.30
8 LCMVV (cm)
6.66
5.45
5.87
7.17
8.92
8.46
104
The above results obtained for the right lung source indicate that the source location
estimates are signiﬁcantly erred for both the pressure and velocity signals. Figure 514
shows that neither the direction of the beamshape nor the range estimates are accurate.
Among the three algorithms, LCMV provides the best estimates, although not as accurate
as desired. Also, the results appear to be marginally better for 7] = 1 as compared to 7] = 0.
5.5.4 Discussion
Comparing the source localization algoritlnm performance for the mitral valve (and
right lung) source case to the homogeneous medium case, it is clear that the heterogeneity
and complexity of the thoracic medium in the former is not handled appropriately by the
nearﬁeld propagation model. Also, the sensor locations used for the two cases were
different, which might also contribute to the difference in the performance of the
algorithms.
In order to compare their performance for the same sensor locations, the algorithms
were tested on the homogeneous thoracic medium with sensor locations as illustrated in
Figure 516, similar to that used in the mitral valve source study.
105
Homogeneous Thorax
=1605 mls, c5 = 0 mls,p=1055 kglm3
z axis (m)
0.2 0.1 0 0.1 0.2
x axis (m)
Figure 516 Location of the true source in the homogeneous thoracic medium is indicated
by the white cross. The white circles correspond to the sensor locations used in the source
localization codes.
Figure 517 shows the objective function images obtained for each of the source
localization algorithms using a loss factor n of O, for only the pressure signals at an SNR
of 5 dB. The images in the ﬁgure also indicate the true source location by a white cross
and the estimated location by a white circle. Further, the title of each image provides the
localization error 8 of the corresponding algorithm.
106
MUSICP 2 a =14.1156 cm
0.05
0
0.05
0.1
0.1 0 0.1
CFBP : e = 5.5036 cm
'0.1 o 0.1
Figure 517 Objective function images of all three source localization algorithms using 11
= 0 obtained on only the pressure input signals at an SNR of 5 dB irn a homogeneous
thoracic medium. Both x and z axes of all images are in m. White cross indicates the true
source location while the white circle corresponds to the estimated source location.
The localization error 8 in cm of each of the algorithms using I] = 0 is tabulated
against SNR in Table 56 for the pressure input signals.
107
Table 56 Localization error 8 in cm versus SNR in dB obtained in all three source
localization algoritlnms using I] = 0, for only pressure input signals in a homogeneous
thoracic medium.
SNR (dB) 15 10 5 0 5 10
r, 0 0 0 0 0 0
aMUSICP(cm) 14.12 14.12 14.12 14.12 14.12 14.12
eCFBP(cm) 6.08 6.08 5.50 13.77 8.39 13.94
eLCMVP(cm) 14.33 14.33 14.12 14.12 14.12 14.12
The above results, at ﬁrst glance, seem quite unexpected. Comparing the source
location estimates for pressure signals in Tables 51 and 56, we notice that for the same
homogeneous medium, a change of sensor locations has resulted in signiﬁcant
performance degadation. Although the directions of the beamshapes seem to be fairly
accurate in Figure 517, the range estimates are highly inaccurate.
This leads us to conclude that the nearﬁeld propagation model as deﬁned in section
(5.1.2) does not accurately portray the signals obtained using the FDTD codes, even for
the homogeneous case. One primary reason for this discrepancy can be attributed to the
source model. The nearﬁeld propagation model is in theory applicable to only a radial
omnidirectional point source. However, in all the F DTD simulations, the source was
assumed to be an f; directional point source.
Thus, in conclusion, two major factors that seem to impair the nearﬁeld propagation
model based source localization algorithms are:
(i) heterogeneity and viscous nature of the medium
(ii) nonradial, non omnidirectional source orientation
108
One of the ways to counter these two issues is to employ a true viscoelastic model of the
thorax obtained via the sourcetosensor transfer functions described in Chapters 3 and 4
in sections (3.7.2) and (4.3). This approach is undertaken in the following section.
109
5.6 Transfer Function based Source Localization Algorithms
As mentioned in the concluding remarks of the previous section, a transfer ﬁmction
based model is likely to provide better results as it inherently incorporates the effects of
complex inhomogeneities and viscosity of the various tissues into its model. Furtlner, the
transfer functions are obtained from simulations that employ the same source model as
used in generating the sensor signals.
Four different transfer function based source localization algorithms are considered.
They are brieﬂy described in the subsections below. Each of algorithms employ the
sourcetosensor transfer functions, say hk1(a))=[x1(a)) / ﬁk(0))] in their algorithms,
where k and 1 correspond to the source and sensor location indices, respectively, (0
corresponds to the angular frequency, x1(a)) is the Fourier spectrum of the 1th sensor
signal and fzk( (0) is the Fourier spectrum of the kth fz source function. For every source
location k, the viscoelastic FDTD simulations provide the transfer functions hk1(a)) for
all the sourcetosensor combinations. In all 56 such sinnulations corresponding to 56
distinct source locations are carried out, each simulation employing N = 19 sensor sites,
thus providing 56 x 19 transfer ﬁnnctions hk1(0)), k = 1,2,...56 and l = 1,2,. . .N. Figure 5
18 shows the locations of all the sources and sensors considered.
110
Source and Sensor Locations
2 axis (m)
O.15 0.1 0.05 0 0.05 0.1 0.15
xaxis (m)
Figure 518 All source and sensor locations considered for the transfer functions
5.6.1 Transfer Function  Source Waveform Energy (SWE)
The name “Source Waveform Energy (SWE)” stems from the fact that this algorithm
attempts to ﬁnd the location that maximizes the energy of the source waveform estimate.
The source function energy estimate for one particular location k is obtained as follows:
(i) The source function fzk1(a)) is computed using the inverse of the transfer
function hk1(a)) for every sensor location I as follows:
111
(ii)
(iii)
the ind
5.6.2
algon
VeCto
2.114) = :53) w
(61)
(ii) It is then ﬁltered using a bandpass ﬁlter with cutoff ﬁ'equencies of 20 Hz and
700 Hz and transformed into time domain using the inverse Fourier transform
to yield fzkl (t).
(iii) Each of the 1 time domain source function estimates are tlnen normalized.
Finally, the median signal from the 1 time domain signals is computed and its
energy Pkis stored as the source waveform energy for the location k.
x 2
median kal(t) "
P = x
k I max fzkl (t )"
(62)
The above procedure is repeated for all k. The source location is ﬁnally estimated as
the index k that maximizes Pk.
5.6.2 Transfer Function  MUSIC
The transfer function based MUSIC algorithm is derived from the original MUSIC
algorithm described in section (5.2.2). The only difference is that the nearﬁeld steering
vector a(r, 6, Cl) is replaced by the transfer function hk(a)) = [hk1(a)), hk2(a)),
...h/WM)”T in objective function equation (54). Since the transfer function is in the
frequency domain, the eigen vectors corresponding to the signal and noise subspace are
112
.rr’
also computed in the frequency domain [73]. The objective function of the transfer
ﬁmction based MUSIC thus becomes:
Pk(w)= 1 2
241141414
i=K +1
(63)
where vi(a)), i = K+1, K+2,. . .,N correspond to the (NK) noise eigen vectors of the data
covariance matrix computed in the ﬁequency domain. The objective ﬁnnction Pk is
computed only at the central ﬁ'equency in the current work. However, a wideband
extension involving all frequencies is trivial.
5.6.3 Transfer Function  CFB
The transfer function based CFB, popularly known as Matched Field Processing
follows the same idea as explained above [22], [73], i.e. the nearﬁeld steering vector
a(r, 6, in) is replaced by the transfer function hk(a))=[hk1(w),hk2(a)), ...hmayf in
objective function equation (56). The objective function in this case is thus given by:
Pk (w) = th (w)R(w)h(w)
(64)
where R(a)) corresponds to the data covariance matrix in the frequency domain.
5.6.4 Transfer Function — LCMV
The same ideology applies to the transfer function based LCMV algorithm too. The
objective function in this case becomes:
113
114
(65)
5.7 Results using Transfer Function based Techniques
This section presents the results of the SWE, MUSIC, CFB and LCMV algoritlnms
employing the transfer function model on the simulated signals for three different source
locations in the thorax. The testing protocol is the same as explained in section (5.4)
except that in this case there are no loss factor I] and velocity 6 parameters to be
optimized. The search for the intrathoracic source location is carried out in the ﬁxed
rectangular xz search space of 0,, X D2 = 22.4 cm X 19.2 cm inside the thorax as
illustrated in Figure 53. However the spatial gidsize in this case is 3 cm X 3 cm.
5.7.1 Mitral Valve Source
A point source located at the mitral valve in the thorax is considered here as shown in
Figure 510 in section (5.5.2). Figure 519 shows the objective function images obtained
for each of the source localization algorithms for both pressure and velocity input signals
at an SNR of 5 dB. The images in the ﬁgure also indicate the true source location by a
white cross and the estimated location by a white circle. Further, the title of each image
provides the localization error 8 of the corresponding algorithm.
115
SWEP : s = 3.6125 cm SWEV : s = 0.67082 cm
0.1 01
0 0
0.1 0.1
0.1 0 0.1 0.1 O 0.1
MUSICP 2 e = 3.6125 cm MUSICV : e = 3.6125 cm
O.1 0.1
0
. 0.1
0.1 0 0.1 0.1 0 0.1
CFBP : a = 4.0804 cm CFBV : a = 4.0804 cm
0.1 0.1
O
O
X X
0 0
0.1 0.1
01 0 0.1 0.1 0 0.1
LCMVP : s = 3.6125 cm LCMVV : s = 3.6125 cm
0.1 01
Figure 519 Objective function images of all four source localization algorithms obtained
on both pressure and velocity input signals at an SNR of 5 dB for a mitral valve source.
Both x and z axes of all images are in m. White cross indicates the true source location
while the white circle corresponds to the estimated source location.
116
The localization error 8 of each of the algorithms is plotted against SNR in Figure 5
20 for both the pressure and velocity input signals.
P V
D SWEP El SWEV
0,25 1* MUSICP4 0,25» 1* MUSICV.
0 CFBP o CFBV
0.2 4 LCMVP 4 0,2 A LCMVV n
g, 015 g, 015
to to i
0.14  0.1 1 i
El Cl
0.05 1 0.05. 
Q B i E i 51 2 S a Q Q Q
Dr D  0. El El El E1 ,
10 0 10 10 0 10
SNR SNR
Figure 520 Plot of localization error 8 in m versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a mitral valve
source.
Table 57 presents the same results as illustrated in Figure 520 in a tabulated form.
117
Table 57 Localization error a in cm versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a mitral valve
SNR (dB) 15 10 5 o 5 1o
aSWEP (cm) 0.67 3.61 3.61 3.61 2.77 6.13
aMUSICP(cm) 3.61 3.61 3.61 3.61 3.61 3.61
aCFBP (cm) 4.08 4.08 4.08 4.08 4.08 4.08
aLCMVP(cm) 3.61 3.61 3.61 3.61 3.61 3.61
aSWEV (cm) 0.67 0.67 0.67 2.77 0.67 6.18
eMUSlCV(cm) 3.61 3.61 3.61 3.61 2.42 2.42
eCFBV(cm) 4.08 4.08 4.08 4.08 4.08 4.08
aLCMVV(cm) 3.61 3.61 3.61 3.61 2.42 2.42
SOUICC.
The above results indicate that the transfer function model employed by the
algorithms result in signiﬁcantly better estimates of the source locations for both the
pressure and velocity signals as compared to the nearﬁeld model based techniques. From
Table 57, it appears that the SWE algorithm provides the best source location estimates.
5.7.2 Right Lung Source
A point source located in the right lung of the thorax is considered here as shown in
Figure 513 in section (5.5.3). Figure 521 shows the objective function images obtained
for each of the source localization algorithms for both pressure and velocity input signals
at an SNR of 10 dB. The images in the ﬁgure also indicate the true source location by a
white cross and the estimated location by a white circle. Further, the title of each image
provides the localization error 8 of the corresponding algorithm.
118
SWEP I 8 = 8.6313 cm SWEV I a = 0.22361 cm
0.1 0.1
0 0
0.1 O ‘ 0.1
0.1 0 0.1 —0.1 0 0.1
MUSICP I e = 0.22361 cm MUSICV I e = 0.22361 cm
0.1 0.1
0 0
0.1 0.1
0.1 0 0.1 —0.1 0 0.1
CFBP : e = 6.1033 cm CFBV : a = 6.1033 cm
O.1
0.1
0 0
O
0.1 0.1
0.1 0 0.1 —0.1 0 0.1
LCMVP I a = 0.22361 cm LCMVV I a = 0.22361 cm
0.1 ‘ 0.1
0 0
0.1 0.1
0.1 0 0.1 0.1 0 0.1
Figure 521 Objective function images of all four source localization algorithms obtained
on both pressure and velocity input signals at an SNR of 10 dB for a right lung source.
Both x and z axes of all images are in in. White cross indicates the true source location
while the white circle corresponds to the estimated source location.
119
Th
23 for
Table
8 (m)
FigUI
The localization error 8 of each of the algorithms is plotted against SNR in Figure 5
22 for both the pressure and velocity input signals and the same results are tabulated in
Table 58
P V
:1 SWEP D SWEV
025 . * MUSICP , 0.25 _ O MUSICV .
0 CFBP o CFBV ..
g, 0.15 g 0.15 8  L
0.14 “  0.1 4
:1 Cl [:1
0.05% 6 o o o O o 0.05 o o o o
E]
0. l 4 A at . 0.4 l A I I 
10 0 10 10 0 10
SNR SNR
Figure 522 Plot of localization error 8 in m versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a right lung
source.
120
algoril
DTCSSL
Table
identi
veloc
5.7.3
the It“
local
The;
Table 58 Localization error 8 in cm versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a right lung
s
source.
SNR (dB) 15 10 5 0 5 10
aSWEP (cm) 8.63 8.63 3.20 8.63 3.11 6.94
aMUSICP (cm) 0.22 0.22 0.22 10.87 12.80 12.80
aCFBP(cm) 6.10 6.10 6.10 6.10 6.10 6.10
sLCMVP (cm) 0.22 0.22 0.22 10.87 12.80 12.80
eSWEV (cm) 0.22 0.22 12.10 12.10 12.49 15.38
eMUSICV(cm) 0.22 0.22 0.22 12.80 12.80 12.80
sCFBV(cm) 6.10 6.10 6.10 13.46 6.10 16.36'
sLCMVV(cm) 0.22 0.22 0.22 12.80 12.80 12.80
The above results indicate that the transfer function model employed by the
algorithms result in signiﬁcantly better estimates of the source locations for both the
pressure and velocity signals as compared to the nearﬁeld model based techniques. From
Table 58, it appears that the MUSIC and LCMV algorithms provide the best (and
identical) source location estimates. Also, the SWE algorithms seem to fare better with
velocity signals as compared to pressure signals for high SNR values.
5.7.3 Left Lung Source
A point source located in the left lung of the thorax is considered. Figure 523 shows
the locations of the true source and sensors used in the source localization codes.
Figure 524 shows the objective function images obtained for each of the source
localization algorithms for both pressure and velocity input signals at an SNR of 10 dB.
The images in the ﬁgure also indicate the true source location by a white cross and the
121
estimated location by a white circle. Further, the title of each image provides the
localization error 8 of the corresponding algorithm.
The localization error 8 of each of the algorithms is plotted against SNR in Figure 5
25 for both the pressure and velocity input signals and the same results are tabulated in
Table59.
Left Lung Source
01 ' A
0.05
E
.2 0
X
(t!
N
0.05
0.1, ..
0.15 0.1 0.05 0 0.05 0.1 0.15
x axis (m)
Figure 523 Location of the true source in the left lung of the thorax is indicated by the
white cross. The white circles correspond to the sensor locations used in the source
localization codes.
122
SWEP : e = 3.4234 cm SWEV : e = 0.5 cm
0.15 ' . . . .
MUSICP : s = 9.5084 cm MUSICV Z a = 9.5084 cm
0.1
0
0.1
0.1 0 0.17 ' o.1 o 0.1
CFBP : a = 3.4132 cm CFBV 1 s = 3.4132 cm
' o.1 o 0.1 ’ o.1 0 0.1
' 0.1 o 0.1
Figure 524 Objective function images of all four source localization algorithms obtained
on both pressure and velocity input signals at an SNR of 10 dB for a left lung source.
Both x and z axes of all images are in m. White cross indicates the true source location
while the white circle corresponds to the estimated source location.
123
P V
0.3 . . 0.3 . .
D SWEP D SWEV
* MUSICP it MUSICV
0.25  0.2 ' r
o CFBP 5 0 CFBV
A LCMVP A LCMV—V
0.2. ~ 0.2* 4
£06» §0w~ «
(O (D
0.1"» a A a A 41‘ 01'4. 4» a a A 41*
1:1 1:1 1:1 :1
0.05 0.051
O O U 0 U U 0 O O O O
0 0, 1:1 1:1 1:1 1:1
10 0 10 1o 0 1o
SNR SNR
Figure 525 Plot of localization error 8 in m versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a left lung source.
Table 59 Localization error 8 in cm versus SNR in dB obtained in all four source
localization algorithms for both pressure and velocity input signals for a left lung source.
SNR(dB) 15 10 5 0 5 10
eSWEP (cm) 3.42 3.42 6.41 3.42 6.41 6.41
aMUSICP(cm) 9.51 9.51 9.51 9.51 9.51 9.51
8CFB—P(cm) 3.41 3.41 3.41 3.41 3.41 3.41
8LCMVP(cm) 9.51 9.51 9.51 9.51 9.51 9.51
sSWEV (cm) 0.50 0.50 0.50 0.50 0.50 6.55
8MUSICV (cm) 9.51 9.51 9.51 9.51 9.51 9.51
eCFBV(cm) 3.41 3.41 3.41 3.41 3.41 3.41
8LCMVV(cm) 9.51 9.51 9.51 9.51 9.51 9.51
124
The above results indicate that the transfer function model employed by the SWE and
CFB algorithms result in good estimates of the source locations for both the pressure and
velocity signals as compared to the MUSIC and LCMV algorithms. From Figure 524, it
appears that the MUSIC and LCMV algorithms although provide inaccurate source
location estimates, their objective function images show a second peak at the true source
location estimate. Also, the SWE algorithms seem to fare better with velocity signals as
compared to pressure signals.
5.7.4 Discussion
Overall, the transfer function based algorithms outperform the nearﬁeld based
methods in all cases. This was an expected result as indicated in the beginning of section
(5.6). We also notice that among the four localization algorithms considered for the
transfer function based approach, the SWE algorithm using velocity signals appear to be
the most consistent and accurate.
The objective function images of LCMV and MUSIC appeared to have a distinct
peak in the region close to the mitral valve for all cases. In some cases, this peak was
stronger than the true source peak thereby leading to an erroneous source location
estimate. The additional false peak around the mitral valve region may be attributed to
strong reﬂections coming off from that region that corrupt the contributions of the
original signal waveforms from the true source. Since the LCMV and MUSIC techniques
rely predominantly on the stronger eigenvectors of the data covariance matrix, it is
possible that the undesirable reﬂections contribute signiﬁcantly towards the signal
subspace, thereby impairing the source location estimates of the algorithms. The SWE
algorithm, however, by virtue of the median operation in its source waveform estimate
125
and due to the fact that it gives equal weightage to the information gained from all
sensors, is not affected as long as there are enough sensors employed in the algorithm.
126
5.8 Methods to improve NearField Propagation Model based
Techniques
As mentioned in the concluding sections of section (5 .5), the non omnidirectional
source orientation could be one of the potential reasons for impairing performance of the
nearﬁeld propagation model based source localization algorithms. In order to tackle this
problem, this section proposes a method to incorporate an approximate ‘fz directional
source” behavior into the nearﬁeld propagation model. Preliminary results supporting
this proposed scheme are also provided.
5.8.1 Distributed Source Model
The idea of employing a distributed source model for beamforming applications has
been studied in the past [74], [75]. The ability to model any source orientation using a
distributed set of point omnidirectional sources can be exploited to yield better source
localization results. In the current work, we propose to model the ﬁ, directional source
using a vertical strip of a few, say Q, closely spaced omnidirectional point sources as
shown in Figure 526. Thus, the nearﬁeld propagation model is now governed by
steering vectors associated with each of these Q omnidirectional point sources.
127
sensor i reference
. 0 sensor .
\. .
x. \
\ . r
rik 0‘ \. 0k
'x. \ p
\ O
. \
\‘x 1
\X :
X i
x t
sourcek
Figure 526 Distributed source model
The above scheme is incorporated into each of the nearﬁeld propagation model
algorithms. Preliminary results are presented in the subsection below.
5.8.2 Preliminary results using Distributed Source Model Techniques
The MUSIC, CFB and LCMV algorithms used the distributed source model as
explained above assuming Q = 5. Results obtained on the homogeneous thoracic medium
and for a mitral valve source are presented below.
5.8.2.1 Homogeneous Thorax
In order to compare the performance of the proposed scheme, the algorithms were
tested on the same homogeneous thoracic medium with sensor locations as illustrated in
Figure 516.
128
Figure 527 shows the objective function images obtained for each of the distributed
source based localization algorithms using a loss factor 17 of O, for both pressure and
velocity input signals at an SNR of 5 dB. The images in the ﬁgure also indicate the true
source location by a white cross and the estimated location by a black diamond. Further,
the title of each image provides the localization error 8 of the corresponding algorithm.
The localization error 8 of each of the algorithms is plotted against SNR in Figure 5
28 for both the pressure and velocity input signals and the same results are tabulated in
Table 510.
129
MUSICP : e = 2.4759 cm MUSICV 2 e = 1.7117 cm
0.05
O
0.05
. 0.1
0.1 0 0.1 O.1 0 0.1
CFBP : e = 2.4759 cm CFBV I a = 1.8028 cm
0.05
0
0.05
. 0.1
0.1 0 0.1 0.1 0 0.1
LCMVP : e = 13.7364 cm LCMVV : e = 14.8772 cm
0.05
0
0.05
0.1
0.1 0 0.1
Figure 527 Objective function images of all three distributed source based algorithms
using I] = 0 obtained on both pressure and velocity input signals at an SNR of 5 dB in a
homogeneous thoracic medium. Both x and z axes of all images are in m. White cross
indicates the true source location while the black diamond corresponds to the estimated
source location.
130
P2n=0 V2n=0
«b MUSIC 1. MUSIC
0.25 o CFB  0.25 o CFB
A A
0.2 LCMV 0.2 LCMV q
15,015 A 0 A A A §o15 A A A A A A
to to
01 0.1
A
0.05 0.05
0 5 41 I "" X 0 9 Ilrl D ID 0
10 0 1‘0 10 0 10
SNR SNR
Figure 5—28 Plot of localization error 8 in m versus SNR in dB obtained in all distributed
source based algorithms using I] = O, for both pressure and velocity input signals in a
homogeneous thoracic medium.
Table 510 Localization error 8 in cm versus SNR in dB obtained in all distributed source
based algorithms using 17 = 0, for both pressure and velocity input signals in a
homogeneous thoracic medium.
SNR (dB) 15 10 5 0 5 10
p, 0 0 0 0 0 0
sMUSICP (cm) 2.48 2.48 2.48 1.97 1.71 2.48
a CFBP (cm) 2.48 2.48 2.48 1.97 15.38 1.08
eLCMVP (cm) 1.01 13.74 13.74 13.74 7.52 13.74
a MUSICV (cm) 1.71 1.71 1.71 1.71 1.71 1.71
e CFBV (cm) 1.80 1.80 1.80 1.71 1.71 3.00
eLCMVV(cm) 13.82 13.74 14.88 14.33 14.88 13.82
131
The above results in Table 510 when compared to those obtained by the nearﬁeld
singlesource model in Table 56 shows a dramatic improvement, especially for the
MUSIC and CFB algorithms.
5.8.2.2 Mitral Valve Source
The point source located at the mitral valve location as discussed in section (5.5.2) is
considered here. Figure 529 shows the objective function images obtained for each of
the distributed source based localization algorithms using a loss factor I] of 0, for both
pressure and velocity input signals at an SNR of 5 dB. The images in the ﬁgure also
indicate the true source location by a white cross and the estimated location by a black
diamond. Further, the title of each image provides the localization error 8 of the
corresponding algorithm.
The localization error 8 of each of the algorithms is plotted against SNR in Figure 5
30 for both the pressure and velocity input signals and the same results are tabulated in
Tables 511 and 512.
132
MUSICP : s = 2.4759 cm MUSICV : e = 3.0806 cm
0.05
0
0.05
0.05
0
0.05
0.1 0.1
0.1 0.05 0 0.05 0.1 0.1 0.05 0 0.05 0.1
CFBP I e = 2.4759 cm CFBV 2 1»: = 3.0806 cm
o.05 x o.05
0 0
0.05 0.05
. 0.1
0.1 0.05 0 0.05 0.1 0.1 0.05 0 0.05 0.1
LCMVP : e = 14.8772 cm LCMVV : e = 13.9445 cm
0.05
0
0.05
0.1
0.1 0.05 0 0.05 0.1
Figure 529 Objective function images of all distributed source based algorithms using I]
= 0 obtained on both pressure and velocity input signals at an SNR of 5 dB for a mitral
valve source. Both x and z axes of all images are in m. White cross indicates the true
source location while the black diamond corresponds to the estimated source location.
133
Pzn=0 V:n=0
0.3 . . 0.3 . .
1» MUSIC 1» MUSIC
025’ o CFB ‘ 025' o CFB ‘
0.2, A LCMV . 0.2. A LCMV ,
£0.15'AAAA AlgonsAAAAAA
to A to
0.11 0.1
005» ' 0  005»
 e ' . ' ' o e ' '
e a o
0» . . or . . i +
10 0 10 10 o 10
SNR SNR
P2n=1 V311=1
4 MUSIC t MUSIC
0.251 o CFB ‘ 025’ o CFB
0.2 A LCMV . 02, A LCMV .
£0.15»; A A A A A‘ £0.15'A A A A A A‘
(ID to
0.1 0.1
O t
0051*.0 63 0.05’°,* or
. o 0 .
o_ . . . or . 1 .
10 0 10 10 0 10
SNR SNR
Figure 530 Plot of localization error 8 in m versus SNR in dB obtained in all distributed
source based algorithms using I] = 0 (top row) and I] = 1 (bottom row), for both pressure
and velocity input signals for a mitral valve source.
134
Table 511 Localization error 8 in cm versus SNR in dB obtained in all distributed source
based algorithms using I] = O, for both pressure and velocity input signals for a mitral
valve source.
SNR (dB)
15
10
5
'l
O
0
0
a MUSICP (cm)
4.30
5.59
2.48
6.32
4.53
e CFBP (cm)
4.30
5.59
2.48
6.32
4.53
s LCMVP (cm)
14.88
12.06
14.88
13.74
14.12
e MUSICV (cm)
3.61
3.61
3.08
1.71
1.71
a CFBV (cm)
3.61
3.61
3.08
1.71
1.08
14.33
13.94
13.94
14.33
13.82
s LCMVV (cm)
‘ Table 512 Localization error 8 in cm versus SNR in dB obtained in all distributed source
based algorithms using I] = 1, for both pressure and velocity input signals for a mitral
valve source.
SNR (dB)
15
10
5
10
'I
1
1
1
1
e MUSICP (cm)
6.08
5.88
3.08
7.43
3.76
5.88
8 CFBP (cm)
4.20
5.28
3.08
5.80
3.76
7.77
8 LCMVP (cm)
14.33
14.33
14.33
14.12
14.12
14.33
a MUSICV (cm)
4.74
4.74
2.51
6.08
4.05
6.08
e CFBV (cm)
4.74
4.74
2.51
2.57
2.01
5.79
8 LCMVV (cm)
14.12
14.12
14.33
14.33
14.12
14.33
135
The results obtained for the distributed source based model in Table 511 and Table
512 when compared to those obtained using the nearﬁeld singlesource model in Tables
52 and 53 appear to have improved signiﬁcantly for the MUSIC and CFB algorithms.
136
5.9 Conclusion
From all the source localization studies conducted in the above sections, it is clear
that modelbased techniques can yield improved results when the model matches the true
conditions. The nearﬁeld propagation model, being based on an omnidirectional point
source in a homogeneous medium, does not perform well to localize the intrathoracic
acoustic sources that were generated using an f; directional source. The transferfunction
based approach, being based on the true viscoelastic forward model, is able to provide
very accurate results. However, ﬁrrther studies need to be conducted to assess their
performance in cases of modelmismatch. The distributed source model approach
provides a scheme to improve the performance of the nearﬁeld propagation model and
may be used in absence of a transfer function model.
137
CHAPTER 6
SUMMARY AND CONCLUSIONS
6.1 Summary
A 2D ﬁnite difference time domain viscoelastic model of the human thorax has been
developed that facilitates the study of sound propagation in the thorax. The model
employs the true human anatomical geometry and accurate estimates of the material
properties of all the tissue structures in the thorax. It thus tracks both the shear and
longitudinal waves and all the associated mode conversions that occur at the various
tissues interfaces in the thorax. The FDTD codes are validated for a simple geometry by
comparing the simulation results to the analytical solutions and to results obtained in a
previously published work. Simulations illustrating the sound propagation from the mitral
valve location to the surface of the thorax are presented. The acoustic sourcetoreceiver
transfer functions corresponding to the various auscultation sites on the thorax are brieﬂy
analyzed. Further, the effects of rib structures, shear elasticity of tissues and dimensions
of the thorax on the sound waves are also studied.
The FDTD viscoelastic model has also been used to study the lung sound propagation
in the human thoracic cavity under pneumothorax conditions. By simulating varying
degrees of severity of the disease, the model assists in determining the key frequency
bands that contain the most information to aid in diagnosis. The work thus lends itself for
development of advanced auscultatory techniques for detection of pneumothorax using
noninvasive acoustic sensors.
138
Finally, three different source localization algorithms based on the traditional near
ﬁeld propagation model are studied and implemented. The localizing ability of these
algorithms to detect an intrathoracic source are assessed and compared. In addition, four
methods based on a transferfunction approach are proposed to alleviate the issues of the
nearﬁeld propagation model based techniques. The results of the transfer fImction based
algorithms are compared to those obtained by the nearﬁeld model to demonstrate the
improved localizing ability of the former. Further, one more strategy employing a
distributed source model is proposed to counter one of the issues of the nearﬁeld
propagation model. Results demonstrating the performance improvement are presented.
139
6.2 Conclusions and Future Work
The 2D model in the present state may be used to simulate and study the various
auscultation sounds associated with certain diseased conditions of the heart or lungs, like
pneumothorax, and thereby provide a limited understanding of the following:
o optimal auscultation sites that contain signiﬁcant diagnostic information
o relation of the sensor signals to the extent and location of the damaged tissue or
the diseased state
0 type of multisensor algorithms that might be suitable for localizing the
pathological sites in the thorax
o variations in sensor signals that can be attributed to changes in the torso geometry
Future efforts could be directed towards developing a 3D model of the torso that
would provide a more accurate representation of these thoracic sounds. In addition,
source localization techniques to localize multiple sources in the thoracic cavity and
methods to estimate the source waveforms could also be considered.
The 2D viscoelastic FDTD model needs to be extended to a 3D model for it to be
used as an advanced diagnostic tool. In addition, incorporation of appropriate sound
generation models to account for the various intrathoracic sound sources would add
immense value to the FDTD models, thereby allowing it to be used to design various
array processing and inversion schemes for diagnosis of various thoracic pathologies.
140
REFERENCES
141
REFERENCES
[1] H. Pasterkamp, S. S. Krarnan, and G. R. Wodicka, “Respiratory sounds  Advances
beyond the stethoscope,” Am. J. Respir. Crit. Care Med. 156(3), 974987 (1997).
[2] A. Sarvazyan, “Audibleﬁequency medical diagnostic methods: Past, present and
future,” J. Acoust. Soc. Am. 117, 2586 (2005).
[3] Y. Bahadirlar and H. O. Giilcﬁr, “Cardiac passive acoustic localization:
CARDIOPAL,” Elektrik, Turkish J. of Electrical Engineering & Computer Science 6(3),
243259 (1998).
[4] N. L. Owsley and A. J. Hull, “Beamforrned nearﬁeld imaging of a simulated coronary
artery containing a stenosis,” IEEE Trans. Medical Imaging 17(6), 900909 (1998).
[5] T. J. Royston, X. Zhang, H. A. Mansy, and R. H. Sandler, “Modeling sound
transmission through the pulmonary system and chest with application to diagnosis of a
collapsed lung,” J. Acoust. Soc. Am. 111, 19311946 (2002).
[6] H. Banks and N. Luke, “Modeling of propagating shear waves in biotissue employing
an internal variable approach to dissipation,” CRSC Tech. Report, NCSU (February
2007)
[7] R. Zalter, H. C. Hardy, and A. A. Luisada, “Acoustic transmission characteristics of
the thorax,” J. Appl. Physiol 18(2), 428436 (1963).
[8] J. J. Faber, “Damping of sound on the chest surface,” Circulation Res. 13, 352358
(1963).
[9] L. G. Durand, J. Genest, and R. Guardo, “Modeling of the transfer function of the
heartthorax acoustic system in dogs,” IEEE Trans. Biomed. Eng. BME32(8), 592601
(1985).
[10] J. Verburg, “Transmission of vibrations of the heart to the chestwall,” Adv.
Cardiovasc. Phys. 5(III), 84103 (1983).
[l l] G. R. Wodicka, K. N. Stevens, and D. C. Shannon, “A theoretical model of sound
transmission in the thorax,” Proceedings of the Annual Int. Conf. of the IEEE
Engineering in Engineering in Medicine and Biology Society, Vol. 1, 248249 (1989).
[12] G. Wodicka, K. Stevens, H. Golub, E. Cravalho, and D. Shannon, “A model of
acoustic transmission in the respiratory system,” IEEE Trans. Biomed. Eng. 36, 925934
(1989)
142
[13] G. R. Wodicka, A. Aguirre, P. D. DeFrain, and D. C. Shannon, “Phase delay of
pulmonary acoustic transmission ﬁ'om trachea to chest wall,” IEEE Trans. Biomed. Eng.
39(10), 10531059 (1992).
[14] G. R. Wodicka, P. D. DeFrain, and S. S. Krarnan, “Bilateral asymmetry of
respiratory acoustic transmission,” Med. Biol. Eng. Comput. 32(5), 489494 (1994).
[15] S. Lu, P. Doerschuk, and G. Wodicka, "Parametric phasedelay estimation of sound
transmitted through intact human lung," Med. Biol. Eng. Comput. 33(3), 293298 (1995).
[16] I. V. Vovk, V. T. Grinchenko, and V. N. Oliynik, “Modeling the acoustic properties
of the chest and measuring breath sounds,” Acoust. Phys. 41, 667676 (1995).
[17] I. V. Vovk, V. T. Grinchenko, and V. N. Oliynik, “Measurement of breath noise
from the surface of human thorax,” Proceedings of the IEEE Instrumentation and
Measurement Technology Conference, Vol. 2, 13301335 (1996).
[18] H. Vermarien and E. van Vollenhoven, “The Recording of Heart Vibrations: A
Problem of Vibration Measurement on Soft Tissue,” Med. Biol. Eng. Comput. 22, 168
178 (1984).
[19] A. H. Leung and S. Sehati, “Sound transmission through normal and diseased human
lungs,” Engineering Science and Education Journal 5(1), 2531 (1996)
[20] Q. Grimal, A. Watzky, and S. Naili, “A onedimensional model for the propagation
of transient pressure waves through the lung,” J. Biomechanics 35(8), 10811089 (2002).
[21] T. J. Royston, Y. Yazicioglu, and F. Loth, “Surface response of a viscoelastic
medium to subsurface acoustic sources with application to medical diagnosis,” J. Acoust.
Soc. Am. 113, 11091121 (2003).
[22] M. B. Ozer, S. Acikgoz, T. J. Royston, H. A. Mansy, and R. H. Sandler, “Boundary
element model for simulating sound propagation and source localization within the
lungs,” J. Acoust. Soc. Am. 122(1), 657671 (2007).
[23] C. Narasimhan, R. Ward, K. L. Kruse, M. Guddati, and G. Mahinthakumar, “A high
resolution computer model for sound propagation in the human thorax based on the
Visible Human data set,” Computers in Biology and Medicine 34(2), 177192 (2004).
[24] H. T. Banks and N. S. Luke, “Modeling of propagating shear waves in biotissue
employing an internal variable approach to dissipation,” Communications in
Computational Physics 3(3), 603640 (2008).
[25] Y. C. Fung, Biomechanics: Mechanical Properties of Living Tissues (Springer
Verlag, New York, 1993).
143
[26] M. Cozic, L. G. Durand, and R. Guardo, “Development of a cardiac acoustic
mapping system,” Proceedings of the I 7th International Conference of the IEEE
Engineering in Medicine and Biology Society, Vol. 2, 10291030 (1995).
[27] James V. Candy, “Modelbased signal processing in acoustics: An overview,” J.
Acoust. Soc. Am. 101(5), 3155 (1997).
[28] Y. Bahadirlar and H. O. Gulcur, “Timeﬁequency cardiac passive acoustic
localization,” Proceedings of the 23"! Annual International Conference of the IEEE
Engineering in Medicine and Biology Society, Vol. 2, 18501853 (2001).
[29] M. Kompis and G. R. ' Wodicka, “Spatial information in simultaneous
multimicrophone recordings of thoracic sounds,” J. Acoust. Soc. Am. 99(4), 24782500
(1996)
[30] M. Kompis, H. Pasterkamp, Y. Motai, and G. R. Wodicka, “Spatial representation of
thoracic sounds,” Proceedings of the 20th Annual International Conference of the IEEE
Engineering in Medicine and Biology Society, Vol. 3, 16611664 (1998).
[31] M. Kompis, H. Pasterkamp, and G. R. Wodicka, “Acoustic imaging of the human
chest,” Chest 120, 13091321 (2001).
[32] N. Owsley, “Passive Acoustic Array Processing of Heart Sounds,” First IEEE
Sensor Array and Multichannel Signal Processing Workshop (2000).
[33] A. M. McKee and R. A. Goubran, “Sound Localization in the Human Thorax,”
Proceedings of the IEEE Instrumentation and Measurement Technology Conference,
Vol. 1, 117122 (2005).
[34] A. M. McKee and R. A. Goubran, “Chest sound pickup using a multisensor array,”
Proceedings of IEEE Conference on Sensors, 780783 (2005).
[3 5] John Scmmlow and Ketaki Rahalkar, “Acoustic Detection of Coronary Artery
Disease,” Annual Review of Biomedical Engineering 9, 449469 (2007).
[36] L. E. Kinsler and A. R. Frey, Fundamentals of acoustics (Wiley, New York, 3rd
edition, 1982).
[37] B. A. Auld, Acoustic ﬁelds and elastic waves in solids (John Wiley & Sons, Inc.,
1973)
[3 8] H. L. Oestreicher, “Field and impedance of an oscillating sphere in a viscoelastic
medium with an application to biophysics,” J. Acoust. Soc. Am. 23, 707714 (1951).
[39] A. Hosokawa, “Ultrasonic pulse waves in cancellous bone analyzed by ﬁnite
difference timedomain methods,” Ultrasonics 44(1), e227e23l (2006).
144
[40] Visible Human Project, wwwnhnnihgw/research/visible/visrble hrgnarrhpgrll (date last
viewed 1/9/09).
[41] A Guided Tour of the Visible Human, wwwmgdsci.org/~lvr1_n/VII/tomhtml (date last
viewed 06/03/09), MAD Scientist Network, Washington University Medical School
(1996)
[42] M. Otani, T. Hirahara, S. Shimizu, and S. Adachi, “Numerical simulation of transfer
and attenuation characteristics of softtissue conducted sound originating ﬁ'om vocal
tract,” Applied Acoustics 70(3), 469472 (2009).
[43] T. D. Mast, L. M. Hinkelman, L. A. Metlay, M. J. Orr, and R. C. Waag, “Simulation
of ultrasonic pulse propagation, distortion, and attenuation in the human chest wall,” J.
Acoust. Soc. Am. 106(6), 36653677 (1999).
[44] A. Kernper, C. McNally, E. Kennedy, S. Manoogian, A. Rath, T. Ng, J. Stitzel, E.
Smith, S. Duma, and F . Matsuoka, “Material Properties of Human Rib Cortical Bone
from Dynamic Tension Coupon Testing,” Stapp Car Crash Journal 49, 199230 (2005).
[45] A. Hosokawa, “Simulation of ultrasound propagation through bovine cancellous
bone using elastic and Biot’s ﬁnitedifference timedomain methods,” J. Acoust. Soc.
Am. 118, 17821789 (2005).
[46] T. D. Mast, “Empirical relationships between acoustic parameters in human soft
tissues,” Acoustics Research Letters Online 1, 37—42 (2000).
[47] M. Fatemi, A. Manduca, and J. F. Greenleaf, “Imaging elastic properties of
biological tissues by lowfrequency harmonic vibration,” Proceedings of the IEEE, Vol.
91(10), 15031519 (2003).
[48] X. Zhang, T. J. Royston, H. A. Mansy, and R. H. Sandler, “Radiation impedance of
a ﬁnite circular piston on a viscoelastic halfspace with application to medical diagnosis,”
J. Acoust. Soc. Am. 109, 795802 (2001).
[49] K. Hoyt, T. Kneezel, B. Castaneda, and K. J. Parker, “Quantitative sonoelastography
for the in vivo assessment of skeletal muscle viscoelasticity,” Phys. Med. Biol. 53, 4063
4080 (2008).
[50] E. L. Madsen, H. J. Sathoff, and J. A. Zagzebski, “Ultrasonic shear wave properties
of soft tissues and tissuelike materials,” J. Acoust. Soc. Am. 74(5), 13461355 (1983).
[51] R. Natarajan, J. Williams, and G. Andersson, “Modeling changes in intervertebral
disc mechanics with degeneration,” J. Bone Joint Surgery Am. 88(2), 3640 (2006).
[52] Viscosity, http://hvpertextbool_(.com/physics/matter/viscositv (date last viewed
06/03/09), The Physics Hypertextbook by Glenn Elert, (2008).
145
[53] J. Gennisson, T. Baldeweck, M. Tanter, S. Catheline, M. Fink, L. Sandrin, C.
Cornillon, and B. Querleux, “Assessment of elastic parameters of human skin using
dynamic elastography,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 51(8), 980989
(2004).
[54] M. J ahed and S. J. LaiFook, “Stress wave velocity measured in intact pig lungs with
crossspectral analysis,” J. Appl. Physiol. 76(2), 565571 (1994).
[55] E. L. Mazuchowski and L. E. Thibault, “Biomechanical properties of the human
spinal cord and pia mater,” ASME Key Biscayne, FL. Presented at the Summer
Bioengineering Conference (2003).
[56] J. P. Bérenger, “A perfectly matched layer for absorption of electromagnetic waves,”
J. Comput. Phys. 114, 185200 (1994).
[57] F. H. Drossaert and A. Giannopoulos, “A nonsplit complex ﬁ'equency shifted PML
based on recursive integration for FDTD modeling of elastic waves,” Geophysics 72(2),
T9Tl7 (2007).
[58] F. Collino and C. Tsogka, “Application of the perfectly matched absorbing layer
model to the linear elastodynamic problem in anisotropic heterogeneous media,”
Geophysics 68, 294307 (2001).
[59] Q. H. Liu and J. Tao, “The perfectly matched layer for acoustic waves in absorptive
media,” J. Acoust. Soc. Am. 102(4), 20722082 (1997).
[60] J. Virieux, “PSV wave propagation in heterogeneous media: velocity—stress ﬁnite
difference method,” Geophysics 51, 889901 (1986).
[61] K. S. Yee, “Numerical solution of initial boundary value problems involving
Maxwell's equations in isotropic media,” IEEE Trans. Antennas and Propagation 14, 302
307 (1996).
[62] C. Schroder and W. Scott, Jr., “On the stability of the FDTD algorithm for elastic
media at a material interface,” IEEE Trans. Geosci. Remote Sensing 40(2), 474481
(2002)
[63] G. Eason, J. Fulton, and I. N. Sneddon, “The generation of waves in an inﬁnite
elastic solid by variable body forces,” Phil. Trans. Royal Soc. London (A) 248, 575607
(1956).
[64] S. Ramakrishnan, S. Udpa, and L. Udpa, “A numerical model simulating sound
pr0pagation in human thorax (Periodical style — Accepted for publication)” Proceedings
of IEEE International Symposium on Biomedical Imaging (2009), to be published.
146
[65] S. Ramakrishnan, S. Udpa, and L. Udpa, “A viscoelastic 2D ﬁnitedifference time
domain model of sound propagation in human thorax (Periodical style  Submitted for
publication)” J. Acoust. Soc. Am., submitted for publication.
[66] H. A. Mansy, T. J. Royston, R. A. Balk, and R. H. Sandler, “Pneumothorax
detection using pulmonary acoustic transmission measurements,” Med. Biol. Eng.
Comput. 40, 520525 (2002).
[67] D. Lichenstein, G. Meziere, P. Biderrnan, and A. Gepner, “The 'lung point': an
ultrasound sign speciﬁc to pneumothorax,” Intensive Care Med. 29, 14341440 (2000).
[68] H. A. Mansy, T. J. Royston, R. A. Balk, and R. H. Sandler, “Pneumothorax
detection using computerised analysis of breath sounds,” Med. Biol. Eng. Comput. 40,
526532 (2002).
[69] S. Ramakrishnan, S. Udpa, and L. Udpa, “A numerical model to study auscultation
sounds under pneumothorax conditions (Periodical style — Accepted for publication)”
Proceedings of the 31'" Annual International Conference of IEEE Engineering in
Medicine and Biology Society (2009), to be published.
[70] B. D. Van Veen and K. M. Buckley, “Bearnforrning: A Versatile Approach to
Spatial Filtering,” IEEE ASSP Magazine 5(5), 424 (1988).
[71] I. A. McCowan, D. C. Moore, and S. Sridharan, “Nearﬁeld Adaptive Beamformer
for Robust Speech Recognition,” Digital Signal Processing 12(1), 87106 (2002).
[72] H. Krim and M. Viberg, “Two decades of array signal processing research: the
parametric approach,” IEEE Signal Processing Magazine 13(4), 6794 (1996).
[73] J. G. Berryman, L. Borcea, G. C. Papanicolaou, and C. Tsogka, “Statistically stable
ultrasonic imaging in random media,” J. Acoust. Soc. Am. 112(4), 15091522 (2002).
[74] Y. R. Zheng, R. A. Goubran, and M. ElTanany, “Robust nearﬁeld adaptive
beamforrning with distance discrimination,” IEEE Trans. Speech and Audio Processing
12(5), 478488 (2004).
[75] Y. Jin and B. Friedlander, “Detection of distributed sources using sensor arrays,”
IEEE Trans. Signal Processing 52(6), 15371548(2004).
147
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
IIIIII[III][IIIIIIIIIIIIIIIJIIIIIIIIIIII