. ‘

u '"m‘n‘filfl‘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 fulfillment
of the requirements for the

Doctoral degree in Electrical Engineering

 

 

ML W

Major ProfessorIs Signature
limit 7, WWI
Q fi

Date

 

 

MSU is an Afi‘innative Action/Equal Opportunity Employer

— —.---o-.-._.—r—uu.--o-o-n--n-o-n-I-l-l-c-

 

 

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 fulfillment 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 non-invasively.
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 inverse-source 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-field 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. Specifically, a finite-difference time-domain (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
confidence 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 financially 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 family-like bonds amongst our group.
Finally, I wish to express my heartfelt gratitude to my parents, my sister and my
brother-in—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 Finite-Difference Time-Domain 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 Near-Field Propagation Model ................................................... 77

5.2 Near-Field Propagation Model based Source Localization Algorithms .................. 80
5. 2. 1 Multiple Signal Classification (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 Near-Field 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 Near-Field 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 3-1 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 3-2 Material properties (density, volume compressibility, volume viscosity, shear
elasticity, shear viscosity) of air surrounding the thoracic cavity. .................................... 23

Table 3-3 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 3-4 Material properties (density, volume compressibility, shear elasticity,
resistance coefficient 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 5-1 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 5-2 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 5-3 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 5-4 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 5-5 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 5-6 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 5-7 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 5-8 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 5-9 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 5-10 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 5-11 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 5-12 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 3-1 Photographic cross-sectional slice of the human thorax acquired from the
Visible Human Project database [40], [41] ....................................................................... 20

Figure 3-2 (a) Cross-sectional slice of human thorax partitioned into distinct anatomical
regions, each tissue color-coded with a unique color (b) Corresponding color-code key to
identify the different tissues .............................................................................................. 21

Figure 3-3 Staggered grid layout for the 2D velocity-stress FDTD scheme indicating the
relative spatial and time indices at which the variables are computed .............................. 34

Figure 3-4 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 3-5 j; components of the forcing functions (a) Gaussian first derivative forcing
pulse (b) Sinusoidal forcing pulse ..................................................................................... 43

Figure 3-6 Comparison of vx and V; components obtained via FDTD to their
correspOnding analytical counterparts for the 1 kHz Gaussian first derivative forcing
pulse (all color-bar units are in m/s) .................................................................................. 44

Figure 3-7 Comparison of vJr and v2 components obtained via FDTD to their
corresponding analytical counterparts for the 1 kHz sinusoidal forcing pulse (all color-bar
units are in m/s) ................................................................................................................. 45

Figure 3-8 Comparison of numerical FDTD with analytical results for a 1 kHz Gaussian
first-derivative 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 3-9 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 3-10 FDTD simulation results for the elastic wave propagation problem defined in
Collino et al. [5 8] — Snapshots of absolute velocity over the entire domain caused due to
an explosive source located at the top lefi corner of the domain ...................................... 48

Figure 3-11 Source (white cross) and receiver (white circles and triangles) locations in
the thoracic geometry. (The axes are in units of m). ......................................................... SO

Figure 3-12 A 300Hz Gaussian first derivative source signal j; used in the FDTD
simulations (a) Time waveform of the forcing signal (b) Corresponding Fourier spectrum.
........................................................................................................................................... 52

Figure 3-13 Magnitude of source-to-receiver transfer functions (TF) for a source at the
mitral valve location — TF plots at receiver #s 4, 8, 12, 16 and 20 in the 0-100 Hz
frequency range ................................................................................................................. 53

Figure 3-14 Magnitude of source-to-receiver transfer functions (TF) for a source at the
mitral valve location — TF plots at receiver #3 4, 8, 12, 16 and 20 in the 100-400 Hz
frequency range ................................................................................................................. 54

Figure 3-15 Magnitude of source-to-receiver transfer functions (TF) for a source at the
mitral valve location -— TF image at all receivers in the 0-100 Hz frequency range ......... 55

Figure 3-16 Magnitude of source-to-receiver transfer functions (TF) for a source at the
mitral valve location - TF image at all receivers in the 100-400 Hz frequency range ..... 56

Figure 3-17 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 first-derivative forcing

source propagates fi'om the mitral valve location. (The axes on each of the subplots are in
units of m and the color-bars have units of m/s) ............................................................... 59

Figure 3-18 Normal component of the velocity v,, at receiver #15 due to a 300 Hz
Gaussian first-derivative forcing source at the mitral valve location under three cases —
normal thorax, rib-less thorax and shear-less thorax ......................................................... 61

Figure 3-19 Normal component of the velocity v,. at receiver #12 due to a 300 Hz
Gaussian first-derivative forcing source at the mitral valve location for three different
thoraxsizes(L,XLz)—35cmX25cm,30cmX200mand40cmX3Ocm ................. 63

Figure 4-1 Anatomical simulation of a pneumothoracic air cavity on a 2D cross-sectional
slice of the human thorax ................................................................................................... 66

Figure 4-2 Left column of images show the cross-sectional 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)-to-sensor(x) transfer functions. ................................................................... 69

Figure 4-3 Comparison of the amplitudes of the source-to-sensor transfer functions
corresponding to varying sizes of the pneumothoracic air pocket .................................... 71

Figure 4-4 Comparison of amplitudes of difference of the transfer fimctions of varying
degrees of pneumothoracic severity and that of a normal thorax ...................................... 72

Figure 5-1 Concept of beamforming ................................................................................. 76

xi

Figure 5-2 Near-field propagation model .......................................................................... 77

Figure 5-3 The outer white box corresponds to the x-z search space used in the
algorithms. The inner white box shows the size of the finest cell used in the search. ...... 85

Figure 5-4 Velocity 0 search values .................................................................................. 86

Figure 5-5 Sensor index and locations of the 19 sensor sites used in the source
localization algorithms ...................................................................................................... 87

Figure 5-6 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 5-7 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 5-8 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 5-9 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 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

Figure 5-11 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 5-12 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 5-13 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 5-14 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 5-15 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 5-16 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 5-17 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 5-18 All source and sensor locations considered for the transfer functions ........ 111

Figure 5-19 Objective firnction 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 5-20 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 5-21 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 5-22 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 5-23 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 5-24 Objective fimction 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 5-25 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 5-26 Distributed source model ............................................................................. 128

Figure 5-27 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 5-28 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 5-29 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 5-30 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'fil'flqh

”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 coefficient in the i-:]' direction
forcing functions in the 1 direction
coefficient of volume compressibility
coefficient of shear elasticity

coefficient of volume viscosity

coefficient of shear viscosity

XV

31

“I

«E‘s-Tm“

558..

CID

fo

(x0, 2 0)
0'

B

angular frequency
wave-number

stretching function in the I direction
scaling factor of the stretching fimction

loss factor of the stretching function

stress-rate tensor

strain-rate tensor
stress-rate tensor components
strain-rate tensor components

thickness of PML
scaling order of the PML loss factor profile

theoretical reflection coefficient of PML

longitudinal wave velocity
center fi'equency of forcing function

x-z location coordinates of source
standard deviation of Gaussian pulse
amplitude of forcing pulse

peak magritude at center frequency
time-shift 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

stress-rate cumulative sum components

strain-rate 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

source-to-receiver 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-field propagation model steering vector of km source

near-field propagation model steering matrix

geometric spreading loss factor

exponential loss factor

sensor data covarianCe matrix

th . . .
k ergen vector of data covariance matrix

near-field propagation model steering vector corresponding to the
location (r, 0, CI) in a 3D polar coordinate search space

near—field 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
sigral-to-noise ratio

localization error

th h .
k SOUI'CC-tO-lt sensor transfer firnctron
th .
k source-to-sensor transfer function vector
th .
z component of the k source fimctron

. 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]. Specifically, a multi-sensor acoustic
chest pad could be an extremely low-cost, powerful, light-weight 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 out-patient home monitoring or battlefield situations.

Recent years have seen many studies that attempt to exploit the spatio-temporal
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 fiequencies, 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-field conditions that prevail. The algorithms
would thus seem to benefit if they incorporate a more accurate physics-based 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
finite-difference 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 non-invasive
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 fields 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-finnction 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
heart-thorax 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 fi'om system nonlinearities.

Verburg’s study [10] was one of the most significant works of the 1980s that
attempted to model the heart sound propagation using visco-elastic 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 scientific 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
fi’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 pick-up transducers on the measured surface sounds was
carefirlly analyzed. It was observed that the presence of a transducer inevitably leads to
the conversion of longitudinal wave into surface waves. This can pose significant
challenges when multi-sensor 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 pick-up
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 fi'om the trachea to various sites on the chest wall over a range of frequencies
based on a cross-correlation 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 fiequencies travel at significantly 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 significant developments in the field 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 fi'om the bronchial airways to the skin
surface. In addition, a finite element model simulating a non-axisymmetric 2D geometry
was used to study the same problem. The complex 2D geometry resulted in increased
non-radial 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 significant although the shear
component becomes less significant 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 X-ray 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 non-parenchymal soft tissue regions composed of fat, muscle and
visceral material. Also, the acoustic response field created by an internal acoustic source
was significantly different fiom that predicted by a fi'ee-field 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 finite difference time domain model developed by Narasimhan et al. [23] used
the actual tlnoracic geometry of the thorax obtained fiom the Visible Human Project to
model the sound propagation in the thorax. However, the model assumes a fluid 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 significantly 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 quasi-linear viscoelastic equation proposed by

Fung [25] in a 2D ring-like 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 justifies 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 significant

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 finite-difference 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 multi-sensor 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 model-based signal processing strategies that can provide more meaningful and
useful results.

One such model based approach discussed in [3] uses the Multiple Signal
Classification (MUSIC) algorithm for localization of acoustic sources in the heart by
employing a near-field propagation model. The main drawbacks of this approach are that
it assumes:

o a homogeneous medium of propagation

0 there is no scattering or reflection 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 time-frequency matrix that provides important
clues on the likely source locations in the heart.

Kompis et al. in [29], [30] and [31] employ a similar simplified model assuming
homogeneous wave propagation and a constant damping factor in their least-squares
estimation approach.

Owsley in [4], [32] uses multiple auscultation recordings for localizing the artery-

blockage sites. A near-field focused beamformer is used to perform imaging of the spatial

10

shear wave energy field 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, significant
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 specifically 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 significant 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 reflections 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 sufficient evidence of the transmission
of shear waves in addition to longitudinal waves in the human torso at the lower audible
range fi'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 fiom 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 fluid flow, equation
(1) would include terms with shear and bulk viscosity. Also, tissues show visco-elastic
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 stress-strain (CS-8) 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 stress-strain 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 benefit to combine these two models into
one framework that allows one to use parameters of either model. If we define fi M and

fly 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 stress-strain 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 coefficients in the 1 direction, T m is

the shear stress and 1&2 is the corresponding resistance coefficient in the x-z 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] define the stiffness of the

17

isotropic medium under consideration arnd are termed as coefficients of volume

compressibility and shear elasticity, respectively. ’12 and [12 define the viscosity of the

medium as pararneterized by the Voigt model and are termed as coefficients of volume
viscosity and shear viscosity, respectively. a) is used to represent the angular frequency

of the wave. The resistance coefficients of stress, originating fi'om the Maxwell’s model
of viscosity, also control the viscous absorptive loss in the medium and are related to the

attenuation coefficients of the tissues. Due to the heterogeneity of the medium, the

density p, Lamé constants (k and p) and the resistance coefficients 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 cross-sectional slice at the mid-
thorax level passing through the heart, as shown in Figure 3-1, is chosen as the 2D image
for the model study. The various anatomical features, inclusive of the organ and tissue
boundaries are identified and the image is partitioned accordingly. Figure 3-2 shows this
partitioned color-coded 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 3-4. 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 3-1 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 3-2 (a) Cross-sectional slice of human thorax partitioned into distinct anatomical
regions, each tissue color-coded with a unique color (b) Corresponding color-code key to
identify the different tissues

21

Table 3-1 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) (Pa-s) (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 Pa-s [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 3-2 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) (Pa-s)
Aira 1.24 147k 0.13

 

#1 = #2 = 0 Pa-8 [42]
aReference 42

Table 3-3 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 Pa-s [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 coefficient 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 chest-air interface. There are numerous ways to
incorporate appropriate boundary conditions at this interface. The most common
approach involves imposing a stress-free condition on the torso surface. However, for a
finite-difference approach using a rectangular gid, imposing a stress-flee 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 finite-difference framework seamlessly.

3.4.1 Concept of PML

The basic idea of PML is to create an additional layer of a reflection-less 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

wave-number, respectively. This represents a non-attenuating plane wave. If the space

~ . w
variable x is replaced by a complex stretched version of it, say x = axx = [1+r-i]x ,
a)

tlnen the plane wave solution modifies 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 non-zero 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 reflection-less property of the PML region [56].

Most PML formulations employing this coordinate stretching concept however,
require an artificial splitting of velocity and stress variables in the FDTD implementation
that decreases the computational efficiency of the codes. In order to avoid the
computational burden caused due to splitting, a non-split 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 stretched-coordinate 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 fi'equency domain. The
recursive integation PML technique employs two auxiliary tensors, the stress-rate tensor
S and the strain-rate tensor E through which the non-split set of equations are obtained.

These tensors are defined 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 definitions 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 +2fl2)§xx +i012§zz

(21)
Transforming these back into time-domain, and following the same procedure for all the

equations in (11) — (15), we obtain the following velocity-stress 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 definition of the
complex stretching firnction equation (16) into equations (18) -— (19) and transforming

them into time-domain. 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 non-PML regions. A value of a1 = 1 and w) = O is chosen in the

computational non-PML region, and a positive w; in the PML region to introduce the

loss. The reflection-less 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 reflection cannot be practically achieved. In order to reduce

the numerical reflections in the PML region, the loss factor W] in this region is usually

chosen to have a tapered profile 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 reflection coefficient and CI, is the longitudinal wave velocity in the PML

region [56], [58], [59]. The amount of reflection 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 flow haemodynannics. Similarly, modeling the respiratory breath sounds is
equally or more challenging as it involves complex air flow 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 first 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 fir is set to zero.

P -(t-to)2 1
1 2

J2na

 

fz(x,z,t)=B(t—t0) e 20 5(X-XO)§(Z-ZO) VtZO

 

 

=0, Vt <0
fx(x,z,t) = O, Vt

(30)

where the standard deviation 0' and amplitude B are defined such that the peak magnitude
at the center fiequencyfo in the Fourier domain is kept constant. The forcing function f;

in the fi'equency domain is given by:

30

—a'2a)2 . 7!
N 2 2 flkao +3)
fz(x,z,a))=BO' we e 6(x—x0)6(z-zo)

 

(31)
Thus, by defining 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 time-shift 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)

 

 

 

_ -(t-to)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 define the velocity and stress variables at time t
= 0. Since the forcing function as defined 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 Finite-Difference Time-Domain Implementation

There are a slew of finite-difference time-domain (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 leapfiog 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 time-step while handling the

discontinuities of material properties at the interface boundaries.

3.6.1 Discretization scheme

A staggered gid leap-frog center-differencing scheme [61] was used to discretize the
momentum and constitutive relation equations described in (22) - (26). For this 2D
velocity-stress formulation of the viscoelastic equations, the velocity and stress variables

are used in a staggered gid as shown in Figure 3-3. 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 x-axis and Z-axis

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 3-3 Staggered grid layout for the 2D velocity-stress 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 11-1]
xx 2], 2 xx 2], 2
(ialfljzm +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 second-order accuracy as the other terms

involving the center differencing, and thus does not degrade the overall accuracy.

Similarly, time derivatives of the strain-rate tensor components, i.e. the Voigt’s viscosity

35

terms in equations (24) — (26) are discretized using a three-point stencil to provide
second-order 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>-Tn<w>l

2 Ax

wx . 1 .
—S +—, ,0
2 ”G 2’)

‘A‘ k-l 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 stress-rate and strain-rate components, and by

defining the stress-rate sums and strain-rate sums as follows,

Qxx(i+-:—,j,k-1)= 19,211,,“ [1+],13m J
. 1 . . 1 .
=Q:(ll+§,j,k-2J+Wx5xx[l+-2-,_],k-1)

(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
—Atflxx[i+%,j,k-l)
. 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 — (i-lle-l-J
_ x 2, , 2 ..
-At‘I’xx(i,j,m+%]
‘I’ [i 'm+-3—)-‘P (i 'm+l)+wE [i 'k+—1-)
xx 9.]: 2 fl 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 time-stepping 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 Courant-Levy condition that imposes a limit on the smallest time-step

allowed. According to the Courant condition,

 

 

cpAt\/ 1 2 + l 2 $1
(Ax) (A2)
(43)
must be satisfied 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
ficp,max

 

(44)
Equation (44) is thus used to decide on the time-step for the FDTD scheme.

However, the Courant limit on the smallest time-step need not be sufficient to
guarantee stability in a heterogeneous medium. Schroder [62] shows that the smallest
time-step required for a stable solution in elastodynamic wave propagation problems can
be more restrictive than the Courant condition as defined 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 sufficient for stability in any medium, the material

properties are averaged at the interface regions, especially at boundaries where the

properties differ significantly. 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 simplified 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 firnctions 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 simplified case. We consider a

homogeneous, infinitely 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 3-4 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

time-step 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 3-5(a) shows a first 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 first case while Figure 3-5(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 3-6 and
Figure 3-7 [64]. Figure 3-6 compares the results obtained for the Gaussian first derivative

forcing pulse while Figure 3-7 compares those for the sinusoidal forcing pulse.

Figure 3-8 compares the normalized vz component obtained using the analytical

expression and the FDTD simulation for the 1 kHz Gaussian first 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 reflections are observable in Figure 3-8 [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 3-5]; components of the forcing functions (a) Gaussian first derivative forcing
pulse (b) Sinusoidal forcing pulse

43

Receiver

er#

Receiv

vx-FDID vz-FDTD

. 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 3-6 Comparison of Vx and v2 components obtained via FDTD to their
corresponding analytical counterparts for the 1 kHz Gaussian first derivative forcing
pulse (all color-bar units are in m/s)

vx-FDTD vz-FDTD

   

20
2
E15
3 o
310
0
fl: 5 _2
1 2 3 4 5 1 2 3 4 5
Time (ms) Time (ms)
vx -Analytical 122-Analytical
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 color-bar

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 3-8 Comparison of numerical FDTD with analytical results for a 1 kHz Gaussian
first-derivative point forcingfiz 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 3-9 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 3-9 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 lefi comer of the domain

47

 

t = 16.7ma

Figure 3-10 FDTD simulation results for the elastic wave propagation problem defined 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
cross-sectional 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 time-step was
chosen to be 0.9 times the Courant limit as defined in equation (44). The forcing function

fz was a first derivative of a Gaussian pulse with a center fi'equencyfi) 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 3-11 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 figure.

49

 

Figure 3-11 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 fimctionfz and its frequency spectrum are shown in Figures
3-12(a) and 3-12(b) respectively. The source-to-receiver transfer fiinctions (TF) were
computed for each of the receiver signals. The normal component of the velocity v” was

used for the receiver signals. Figure 3-13 and Figure 3-14 show the magnitude of the

transfer functions of receiver #5 4, 8, 12. 16 and 20 in the frequency range of 0-100 Hz

50

and 100-400 Hz, respectively. Figure 3-15 and Figure 3-16 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 fimction : 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 3-12 A 300Hz Gaussian first 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

|Ve|ocity-Force TF(00)|

A
I

 

 

 

20 4'0 60 80 00
Frequency (Hz)

Figure 3-13 Magnitude of source-to-receiver transfer functions (TF) for a source at the
mitral valve location — TF plots at receiver #5 4, 8, 12, 16 and 20 in the 0-100 Hz
frequency range

53

 

 

—l
1

.0
m

 

 

 

.9
m

p
a.

|Ve|ocity-Force TF(<0)|

0.2 ~

 

 

 

 

100 150 200 250 300 350 400
Frequency (Hz)

Figure 3-14 Magnitude of source-to-receiver transfer functions (TF) for a source at the
mitral valve location — TF plots at receiver #5 4, 8, 12, 16 and 20 in the 100-400 Hz
frequency range

54

|Ve|ocity-Force TF(co)|

   

20
11?
E 40
>.
O
C
d)
g. 60
9
U.

80

100

20 40 60
Receiver#

Figure 3-15 Magnitude of source-to-receiver transfer functions (TF) for a source at the
mitral valve location — TF image at all receivers in the 0-100 Hz frequency range

55

|Ve|ocity-Force TF(o))|

200

250

Frequency (Hz)

300

 

350

 

10 20 30 40 50 60 70
Receiver#

Figure 3-16 Magnitude of source-to-receiver transfer functions (TF) for a source at the
mitral valve location — TF image at all receivers in the 100-400 Hz frequency range

From Figure 3-13 and Figure 3-14, it is clear that the heart~to-chest acoustic transfer
functions at different receiver locations along the front surface of the thorax differ
significantly, although, all of them have most of their energy content lying in the O-80Hz

frequency range. In Figure 3-14 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 fiom the image in Figure 3-16 across all the 20 receivers on the fiont
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 0-40 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 modified 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 cross-sectional slice in the lungs. Since
the frequency range of lung sounds extend well beyond 1 kHz, the center fi'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 3-17 shows the vx and vZ component

57

distributions in the thorax at five time snapshots as the Gaussian first-derivative pulse

source waveform propagates fiom 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 wave-front 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 reflections at the heart-lung 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 ~0-4 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 3-17 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 first-derivative
forcing source propagates from the mitral valve location. (The axes on each of the

subplots are in units of m and the color-bars 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 modifications to the parameters. The
study, whilst aids understanding the effects of such anatomical variations, also justifies

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 3-1 to Table 3-4 that the ribs and bones have distinctly different
material properties as compared to the other tissues. One thus expects their presence to
cause a significant 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

coefficient of shear elasticity ,u] was set to zero in the entire domain to simulate the

shear-less case. The normal component of the velocity V" obtained at receiver #15 for a

normal thorax, rib-less thorax and shear-less thorax are compared in Figure 3-18 [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 significantly,

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 shear-wave contribution in certain tissues like the lung-parenchyma are
insignificant compared to that of the compressional waves, more efficient 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 simplifies 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 3-18 Normal component of the velocity v” at receiver #15 due to a 300 Hz
Gaussian first-derivative forcing source at the mitral valve location under three cases -
normal thorax, rib-less thorax and shear-less thorax

61

3.8.2 Effect of size of thorax

All the simulations in the preceding sections have assumed a fixed 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 fit 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 3-19 [65]. Results indicate that
while the amplitude of the signal corresponding to a larger sized thorax is smaller and
time-delayed compared to that of the smaller sized thorax, the shape of the waveform is
retained atleast upto the first wavefiont, corresponding to the dominant longitudinal wave
peak and the shear wave peak. The secondary reflections 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 3-19 Normal component of the velocity v" at receiver #12 due to a 300 Hz
Gaussian first-derivative 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 lung-collapse
eventually. Since it is a potentially life threatening disease, and is easily treatable when
diagnosed early, there is significant interest in developing a noninvasive, low-cost, quick
and advanced diagnostic scheme. Currently, the gold-standard technique for diagnosis of
pneumothorax is the chest X-ray 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 multi-sensor processing scheme, one could devise a
low-cost, noninvasive and portable solution for accurate diagnosis that may be even used
for out-patient 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 multi-sensor 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 multi-sensor 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 source-to-sensor 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 4-1. 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

defined 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 4-1 Anatomical simulation of a pneumothoracic air cavity on a 2D cross-sectional
slice of the human thorax

66

To model the sound wave propagation in the thoracic cavity, we first discretize the

transverse cross-sectional 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
defined 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 source-to-sensor 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 4-2. 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 4-2 represent the
amplitudes of the source-to-sensor pressure TFs for the corresponding pneumothorax
condition. The image and TF plot in the first 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 figure and subsequent figures.

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 4-2 Lefi column of images show the cross-sectional slice of the thorax for varying
sizes of the pneumothoracic air pocket in the lung. The labels on the lefi side indicate the
size of the air pocket. The right column of plots represent the corresponding amplitude of

the source(o)-to-sensor(x) transfer functions.

69

In order to firrther study and compare these TFs, Figure 4-3 plots these on one
common axis. We notice that at frequencies around 50-100 Hz there isn’t any significant
difference in the amplitudes of the TFs. However, the frequency range of 350-475 Hz
shows a distinct trend of decreasing amplitudes with increasing degree of severity of
pneumothorax. Figure 4-4 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 350-475
Hz fi'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 significantly distant from the air cavity do not show such prominent
deviation. Thus, by combirning information from various sensor sites, one could devise a

multi-sensor 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 4-3 Comparison of the amplitudes of the source-to-sensor 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 test-bed 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 significantly enhanced vvitln a multi-sensor
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 multi-sensor 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 intra-thoracic 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 modified 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
fi'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 5-1. 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 5-1 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
field 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-field problem, which is far more complex. In recent times, there
have been attempts to modify the far-field models to make them suitable to near field
propagation problems [3], [4]. However, the extension from far-field to near-field 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-field 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 Near-Field Propagation Model
The near-field propagation model assumed by the algorithms is based upon a

spherical wave-front 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 5-2.

sensor 1' reference

sensor
0 9 e

 

\O
O ‘. O
s. .
x. \
r,\e \rUk
1k \ .
‘. \
\ O
. \
\..
sourcek

Figure 5-2 Near-field 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 zero-mean 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[_ fl—{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-field propagation model [3], [71]. In equation (52), we notice that

the model contains three explicit terms that are described below:
(i) The first term is related to the geometric spreading loss due to the spherical
wave expanding from an omnidirectional point source. The factor fl, 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 fi'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 defined

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 Near-Field Propagation Model based Source Localization

Algorithms

The three algorithms employing the near-field propagation model described in the
previous subsection considered are:

i) MUltiple SIgnal Classification 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 subspace-based estimation technique that uses

the eigen structure of the sensor data-covariance 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 eigen-vectors (with an eigen-value 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 eigen-vector (i.e. k > K),
tlnen the following hold good:

Rvk = ozvk

H 2 .
=> (ARSA + o I)vk = o’zvk fiom equation (53)
H
=> ARSA V], = 0
=> AHvk = 0 (since Rs and A are botln of rank K)

This implies that the (N — K) noise eigen-vectors are ortlnogonal to the K source

direction vectors, and lie in the null-space of A. This forms the basis for most of the
eigen-vector based algorithms. The MUSIC algorithm uses this concept and estimates the

source locations by identifying the peaks of a function defined as follows:

1

N
vaa(r,t9,¢4

k=K +1

P(r, 6, ¢) =

 

2

 

 

(54)
where a( r, 0, D) is the near-field 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 defined 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 non-uivial solution the weight vector is restricted
to have a norm of 1. Substituting the sensor signal model as defined 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’ ¢>Ra<n 44
a” (r, 9, ¢)a(r, 49, ¢)
(56)

Thus, the source location for the CF B is obtained by finding the (r, 0, C1) in the search

space that maximizes P( r, 0, E1) in equation (56).

5.2.3 Linearly Constrained Minimum Variance Beamformer (LCMV)
The linearly constrained rrninimum variance beamformer, also termed as the Capon’s
beamformer, finds the optimal weight vector W by mirnirrnizing the variance in the output

while subject to constraints that preserve the signals from certain desired directions and

mirnimize contributions from noise or interfering sources arriving from other directions

[71], [72]. If C = [ac 1, dog, . . . acq] are the q constraint steering vector directions and g

= [g 1, g2, ...gq]r are the corresponding constraint values, then the LCMV formulates the

problem as:

82

min WHRW subject to CHW = g

w
(57)
The optimum weights are then provided by:
l
w=R4CkHR4CIg
(53)
The output power for this weight vector thus becomes:
1
P(r,t9,¢) = gH [CHR-IC] g
(59)

One of the most popular approaches of LCMV is to minimize the output power

contributed by noise and any signals coming from directions other than (r, 0, [1) while
maintaining a fixed urnity gain in the (r, 6’, D) direction, i.e. the constraint matrix C is

set to a(r, 0, C1) and g to 1. In such cases, equation (59) becomes:

 

l
P , t9, =
(r ‘1’) a” (r, (9, ¢)R—1a(r, 61,45)

(60)
Thus, the source location for the LCMV is obtained by finding the (r, 0, C1) in the

search space that maximizes P( r, 6, D) in equation (60).

83

5.3 Algorithm Parameters

In the sections to follow, each of the above discussed source localization algorithms
are tested on the simulated thorax signals obtained using the viscoelastic FDTD model in

Chapter 3. The search for the intra-thoracic source location is carried out in a fixed
rectangular x-z search space of Dx X Dz = 22.4 cm X 19.2 cm inside the thorax that is
discretized into tiny rectangular cells governed by the spatial gid-size parameter. Figure
5-3 shows this search domain and the corresponding gid size. The steering vector a(x,
2) corresponding to the center of each of these cells is used to compute the power P(x, 2)

for each (x, 2) cell. For all the simulations, the spatial gid size was set to 8 mm X 8 mm.

A finer gid size could yield more accurate results in cases where the true source location
does not coincide with the center of a cell. However, a search over a finer gid space

requires more computational time.

84

Search Space and Grid Size

   
   

z axis (m)
C

-0.15 -0.1 -0.05 0 0.05 0.1 0.15
xaxis (m)

Figure 5-3 The outer white box corresponds to the x-z search space used in the
algorithms. The inner white box shows the size of the finest cell used in the search.

The three parameters governing the near-field propagation model as described in

section (5.1.2) were optimized as follows:

(i) The geometric spread factor ,8 was set to 2 as this yielded the best results for
all cases.
(ii) The loss factor 17 was searched over two values of 0 and 1. For cases

pertaining to a homogeneous medium, I] was set to 0.

85

(iii) The velocity 6 was searched over a wide range of values fiom 5 m/s to 3500

m/s as plotted in Figure 5-4. For cases, pertaining to a homogeneous medium,

tlne search was conducted over [100 m/s, 1500 nn/s, 1605 m/s].
3500 ' ' r x-
3000
2500
2000

' 1500 x“
10001 xx -

500- xxx -

xxxx

 

x
X
r r

X
r

Veloczty 6 (m/s)
X

 

 

x
0 ~xx>o<”0000000“x

0 10 20 30 40
Velocity Index

 

Figure 5-4 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 fi'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 5-5 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 time-window of the sensor signals used
also dictate the performance of the algorithms. A time-window 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 significantly 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 first derivative pulse as the f},

sourcing function.
White Gaussian noise of varying Signal-to-Noise Ratio (SNR) levels, ranging

from 15 dB to -10 dB, was added to artificially introduce noise into the sensor

signals. Figure 5-6 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 filtered using a 7th order low-pass Chebyshev filter

with the stop-band ripple 20 dB down and with a cut-off fi'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 x-z 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 5-6 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 filter to eliminate any noisy peaks in the function.

(vi) The x-z location of the maximum intensity in the smoothened objective

function for every I] and velocity 0 is stored. The median x-z location across

the range of velocities is used as the final 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 MUSIC-P while that obtained on velocity v" waveforms is titled

as MUSIC-V. Similar convention is adhered to for the remaining two localization

algorithms (CFB and LCMV).

91

 

5.5 Results using the Near-Field Propagation Model based

Techniques

This section presents the results of the MUSIC, CFB and LCMV algorithms
employing the near-field 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 shear-less and viscous-less. In other

words, shear wave contributions are absent and the medium is loss-less. Figure 5-7 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 5-7 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 5-8 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 figure 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

MUSIC-P I a = 0.92195 cm MUSIC-V 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

CFB-P 1 e = 1.7117 cm

-0.05
0
0.05

0.1 .
-0.1 0 0.1 -0.1 0 0.1

LCMV-P I e = 0.92195 cm LCMV-V Z a = 5.5036 cm
9

   

-0.05 x

   

'-0.1 0 0.1 01 0 0.1

Figure 5-8 Objective firnction 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 5-9 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.05-tastes.
A A Q 9 9 O
0' . . . 1 0— . . .
-10 0 10 -10 0 1O
SNR SNR

Figure 5-9 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 5-1 presents the same results as illustrated in Figure 5-9 in a tabulated form.

Table 5-1 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 MUSIC-P (cm)

1.97

1.97

0.92

1.71

15.16

8.19

 

s CFB-P (cm)

2.69

2.69

1.71

4.94

12.59

7.83

 

s LCMV-P (cm)

2.69

2.69

0.92

0.92

1.22

5.59

 

s MUSIC-V (cm)

5.59

5.59

5.59

5.50

8.81

5.59

 

a CFB-V (cm)

5.50

5.50

5.50

5.53

8.81

3.45

 

 

 

a LCMV-V (cm)

5.50

 

5.50

 

5.50

5.59

 

5.50

 

 

5.59

 

 

95

The above results indicate that the near-field model assumed by the algorithms fit the
pressure signals better than the velocity signals. From Figure 5-8, it appears that although
the direction of the beam-shape is fairly accurate, the range estimate obtained for the
velocity signals is erred. For the pressure signals, the source location estimates are
significantly 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 5-10
shows the locations of the true source and sensors used in the source localization codes.

Figure 5-11 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 figure 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

MUSIC-P : a = 14.1156 cm MUSIC-V I e = 11.8004 cm

—0.05
0
0.05

0.1 .
-0.1 O 0.1 -0.1 0 0.1

CFB-P : a = 4.1049 cm CFB-V : e = 5.883 cm

 

    

O

  

0.1 - .,
-0.1 0 0.1 -0.1 0 0.1

LCMV-P I e = 7.7389 cm LCMV-V I e = 11.0368 cm

   

30.1 0 0.1

Figure 5-11 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 5-12 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
$50-15": 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 5-12 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 5-2 and 5-3 present the same results as illustrated in Figure 5-12 in a tabulated

form.

99

 

‘_.:IiI.L

Table 5-2 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 MUSIC-P (cm)

14.33

14.33

14.12

14.12

14.33

14.33

 

a CFB-P (cm)

5.17

4.60

4.05

4.60

5.37

9.54

 

a LCMV-P (cm)

10.82

9.84

14.12

14.12

14.12

14.12

 

a MUSIC-V (cm)

14.12

14.12

14.12

14.12

11.80

14.12

 

e CFB-V (cm)

6.86

6.86

6.86

5.53

0.61

6.86

 

e LCMV-V (cm)

 

 

14.33

 

14.12

 

12.06

 

14.12

 

14.12

 

14.12

 

 

 

Table 5-3 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 MUSIC-P (cm)

14.33

14.12

11.80

11.80

14.33

14.33

 

a CFB-P (cm)

4.60

4.10

3.47

3.72

4.60

9.54

 

a LCMV-P (cm)

9.27

7.74

14.12

14.12

14.12

14.12

 

r: MUSIC-V (cm)

11.80

11.80

11.80

12.57

7.00

14.12

 

e CFB-V (cm)

5.88

5.88

6.58

4.74

2.51

6.58

 

 

e LCMV-V (cm)

 

14.33

 

11.04

 

9.53

 

14.12

 

12.81

 

14.12

 

 

The above results indicate that the near-field model assumed by the algorithms result
in significantly erroneous estimates of the source locations for both the pressure and

velocity signals. From Figure 5-11, it appears that neither the direction of the beam-shape

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 5-13 shows

the locations of the true source and sensors used in the source localization codes.

Right Lung Source

“—5.. ,
5.:Tfi5'5, A" -

 

z axis (m)
0

 

 

I:

-0.15 -0.1 -0.05 0 0.05 0.1 0.15
xaxis (m)

Figure 5-13 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

MUSIC-P 2 e = 6.5765 cm MUSIC-V I e = 7.1701 cm

-0.05
0
0.05

. 0.1
-0.1 0 0.1 -0.1 0 0.1

CFB-P : e = 8.6377 cm CFB-V : 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 5-14 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 5-14 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 figure 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»; . . o-g . .
-10 0 10 -10 0 10
SNR SNR
P1n=1 V2n=1
4 MUSIC 1* MUSIC
0-25' 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 5-15 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 5-15 for botln the pressure and velocity input signals. Tables 5-4

and 5-5 present the same results as illustrated in Figure 5-15 in a tabulated form.

Table 5-4 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 MUSIC-P (cm)

8.46

8.46

8.46

8.46

8.92

7.78

 

e CFB-P (cm)

10.15

10.15

9.95

9.95

10.74

9.44

 

e LCMV-P (cm)

6.61

6.03

6.58

8.46

8.30

7.71

 

c MUSIC-V (cm)

8.46

8.46

8.46

8.46

8.92

8.46

 

e CFB-V (cm)

12.18

12.18

12.18

8.05

12.68

12.30

 

 

a LCMV-V (cm)

 

8.30

 

7.80

 

7.16

 

8.46

 

8.92

 

8.46

 

 

 

Table 5-5 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 MUSIC-P (cm)

6.58

6.58

6.58

8.46

8.92

6.61

 

e CFB-P (cm)

8.39

8.64

9.17

9.17

10.02

8.73

 

s LCMV-P (cm)

5.40

4.84

5.14

8.46

6.23

5.87

 

8 MUSIC-V (cm)

7.17

7.17

7.17

6.58

8.30

8.46

 

a CFB-V (cm)

11.39

11.39

11.39

5.45

12.68

12.30

 

 

8 LCMV-V (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 significantly erred for both the pressure and velocity signals. Figure 5-14
shows that neither the direction of the beam-shape 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-field 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 5-16, 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 5-16 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 5-17 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 figure 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

MUSIC-P 2 a =14.1156 cm

-0.05
0
0.05

0.1
-0.1 0 0.1

CFB-P : e = 5.5036 cm

 

 

'-0.1 o 0.1

 

Figure 5-17 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 5-6 for the pressure input signals.

107

Table 5-6 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

 

aMUSIC-P(cm) 14.12 14.12 14.12 14.12 14.12 14.12
eCFB-P(cm) 6.08 6.08 5.50 13.77 8.39 13.94
eLCMV-P(cm) 14.33 14.33 14.12 14.12 14.12 14.12

 

 

 

 

 

 

 

 

 

 

 

The above results, at first glance, seem quite unexpected. Comparing the source
location estimates for pressure signals in Tables 5-1 and 5-6, we notice that for the same
homogeneous medium, a change of sensor locations has resulted in significant
performance degadation. Although the directions of the beam-shapes seem to be fairly
accurate in Figure 5-17, the range estimates are highly inaccurate.

This leads us to conclude that the near-field propagation model as defined 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-field propagation model is in theory applicable to only a radial

omni-directional 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-field propagation
model based source localization algorithms are:

(i) heterogeneity and viscous nature of the medium

(ii) non-radial, non omni-directional 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 source-to-sensor 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 fimction
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 briefly described in the subsections below. Each of algorithms employ the

source-to-sensor transfer functions, say hk1(a))=[x1(a)) / fik(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 source-to-sensor 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 finnctions 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 5-18 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 find 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 filtered using a band-pass filter with cut-off fi'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 finally 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-field 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

fimction based MUSIC thus becomes:

Pk(w)= 1 2

241141-414

i=K +1

 

 

 

 

(63)
where vi(a)), i = K+1, K+2,. . .,N correspond to the (N-K) noise eigen vectors of the data
covariance matrix computed in the fiequency domain. The objective finnction Pk is

computed only at the central fi'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-field 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 intra-thoracic source location is carried out in the fixed

rectangular x-z search space of 0,, X D2 = 22.4 cm X 19.2 cm inside the thorax as

illustrated in Figure 5-3. However the spatial gid-size 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 5-10 in section (5.5.2). Figure 5-19 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 figure 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

 

SWE-P : s = 3.6125 cm SWE-V : s = 0.67082 cm
-0.1 -0-1

0 0

    

0.1 0.1
-0.1 0 0.1 -0.1 O 0.1
MUSIC-P 2 e = 3.6125 cm MUSIC-V : e = 3.6125 cm
-O.1 -0.1

0

  
  

 

. 0.1
-0.1 0 0.1 -0.1 0 0.1
CFB-P : a = 4.0804 cm CFB-V : 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
LCMV-P : s = 3.6125 cm LCMV-V : s = 3.6125 cm
-0.1 -0-1

   

Figure 5-19 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 SWE-P El SWE-V
0,25- 1* MUSIC-P4 0,25» 1* MUSIC-V.
0 CFB-P o CFB-V
0.2- 4 LCMV-P 4 0,2- A LCMV-V 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 5-20 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 5-7 presents the same results as illustrated in Figure 5-20 in a tabulated form.

117

Table 5-7 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
aSWE-P (cm) 0.67 3.61 3.61 3.61 2.77 6.13
aMUSIC-P(cm) 3.61 3.61 3.61 3.61 3.61 3.61
aCFB-P (cm) 4.08 4.08 4.08 4.08 4.08 4.08
aLCMV-P(cm) 3.61 3.61 3.61 3.61 3.61 3.61
aSWE-V (cm) 0.67 0.67 0.67 2.77 0.67 6.18
eMUSlC-V(cm) 3.61 3.61 3.61 3.61 2.42 2.42
eCFB-V(cm) 4.08 4.08 4.08 4.08 4.08 4.08
aLCMV-V(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 significantly better estimates of the source locations for both the

pressure and velocity signals as compared to the near-field model based techniques. From

Table 5-7, 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 5-13 in section (5.5.3). Figure 5-21 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 figure 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

  
  
  

 

SWE-P I 8 = 8.6313 cm SWE-V 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
MUSIC-P I e = 0.22361 cm MUSIC-V 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

CFB-P : e = 6.1033 cm CFB-V : 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
LCMV-P I a = 0.22361 cm LCMV-V 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 5-21 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 SWE-P D SWE-V
025 . * MUSIC-P , 0.25 _ O MUSIC-V .
0 CFB-P o CFB-V ..
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 5-22 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 5-8 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
aSWE-P (cm) 8.63 8.63 3.20 8.63 3.11 6.94
aMUSIC-P (cm) 0.22 0.22 0.22 10.87 12.80 12.80
aCFB-P(cm) 6.10 6.10 6.10 6.10 6.10 6.10
sLCMV-P (cm) 0.22 0.22 0.22 10.87 12.80 12.80
eSWE-V (cm) 0.22 0.22 12.10 12.10 12.49 15.38
eMUSIC-V(cm) 0.22 0.22 0.22 12.80 12.80 12.80
sCFB-V(cm) 6.10 6.10 6.10 13.46 6.10 16.36'
sLCMV-V(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 significantly better estimates of the source locations for both the
pressure and velocity signals as compared to the near-field model based techniques. From
Table 5-8, 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 5-23 shows
the locations of the true source and sensors used in the source localization codes.

Figure 5-24 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 figure 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

 

 

Table5-9.
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 5-23 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

SWE-P : e = 3.4234 cm SWE-V : e = 0.5 cm

  
    

   

0.15 ' . . . .

MUSIC-P : s = 9.5084 cm MUSIC-V Z a = 9.5084 cm
-0.1

0

 

0.1

 

-0.1 0 0.17 ' -o.1 o 0.1
CFB-P : a = 3.4132 cm CFB-V 1 s = 3.4132 cm

 

      

 

' -o.1 o 0.1 ’ -o.1 0 0.1

   

' -0.1 o 0.1

Figure 5-24 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 SWE-P D SWE-V
* MUSIC-P it MUSIC-V
0.25- - 0.2 ' r
o CFB-P 5 0 CFB-V
A LCMV-P A LCMV—V
0.2. ~ 0.2* 4
£06» §0w~ «
(O (D
0.1"» a A a A 41‘ 0-1'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 5-25 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 5-9 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
eSWE-P (cm) 3.42 3.42 6.41 3.42 6.41 6.41
aMUSIC-P(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
8LCMV-P(cm) 9.51 9.51 9.51 9.51 9.51 9.51
sSWE-V (cm) 0.50 0.50 0.50 0.50 0.50 6.55
8MUSIC-V (cm) 9.51 9.51 9.51 9.51 9.51 9.51
eCFB-V(cm) 3.41 3.41 3.41 3.41 3.41 3.41
8LCMV-V(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 5-24, 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-field 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 reflections 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 eigen-vectors of the data covariance matrix, it is
possible that the undesirable reflections contribute significantly 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 Near-Field Propagation Model based

Techniques

As mentioned in the concluding sections of section (5 .5), the non omni-directional
source orientation could be one of the potential reasons for impairing performance of the

near-field 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-field 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 omni-directional sources can be exploited to yield better source

localization results. In the current work, we propose to model the fi, directional source

using a vertical strip of a few, say Q, closely spaced omnidirectional point sources as
shown in Figure 5-26. Thus, the near-field 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 5-26 Distributed source model
The above scheme is incorporated into each of the near-field 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 5-16.

128

Figure 5-27 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 figure 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 5-10.

129

 

 

 

MUSIC-P : e = 2.4759 cm MUSIC-V 2 e = 1.7117 cm

   

   

 

-0.05
O
0.05
. 0.1
-0.1 0 0.1 -O.1 0 0.1
CFB-P : e = 2.4759 cm CFB-V I a = 1.8028 cm
-0.05
0
0.05
. 0.1
-0.1 0 0.1 -0.1 0 0.1
LCMV-P : e = 13.7364 cm LCMV-V : e = 14.8772 cm
-0.05
0
0.05
0.1
-0.1 0 0.1

 

Figure 5-27 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 5-10 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
sMUSIC-P (cm) 2.48 2.48 2.48 1.97 1.71 2.48
a CFB-P (cm) 2.48 2.48 2.48 1.97 15.38 1.08
eLCMV-P (cm) 1.01 13.74 13.74 13.74 7.52 13.74
a MUSIC-V (cm) 1.71 1.71 1.71 1.71 1.71 1.71
e CFB-V (cm) 1.80 1.80 1.80 1.71 1.71 3.00
eLCMV-V(cm) 13.82 13.74 14.88 14.33 14.88 13.82

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

131

 

The above results in Table 5-10 when compared to those obtained by the near-field
single-source model in Table 5-6 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 5-29 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 figure 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 5-11 and 5-12.

132

 

MUSIC-P : s = 2.4759 cm MUSIC-V : 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
CFB-P I e = 2.4759 cm CFB-V 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
LCMV-P : e = 14.8772 cm LCMV-V : e = 13.9445 cm

-0.05
0
0.05

   

0.1
-0.1 -0.05 0 0.05 0.1

Figure 5-29 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 ‘ 0-25' o CFB ‘
0.2, A LCMV . 0.2. A LCMV ,
£0.15'AAAA Algons-AAAAAA-
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 ‘ 0-25’ 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 5-30 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 5-11 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 MUSIC-P (cm)

4.30

5.59

2.48

6.32

4.53

 

e CFB-P (cm)

4.30

5.59

2.48

6.32

4.53

 

s LCMV-P (cm)

14.88

12.06

14.88

13.74

14.12

 

e MUSIC-V (cm)

3.61

3.61

3.08

1.71

1.71

 

a CFB-V (cm)

3.61

3.61

3.08

1.71

1.08

 

 

14.33

 

13.94

 

13.94

 

14.33

 

13.82

 

 

 

s LCMV-V (cm)

‘ Table 5-12 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 MUSIC-P (cm)

6.08

5.88

3.08

7.43

3.76

5.88

 

8 CFB-P (cm)

4.20

5.28

3.08

5.80

3.76

7.77

 

8 LCMV-P (cm)

14.33

14.33

14.33

14.12

14.12

14.33

 

a MUSIC-V (cm)

4.74

4.74

2.51

6.08

4.05

6.08

 

e CFB-V (cm)

4.74

4.74

2.51

2.57

2.01

5.79

 

 

8 LCMV-V (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 5-11 and Table
5-12 when compared to those obtained using the near-field single-source model in Tables

5-2 and 5-3 appear to have improved significantly 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 model-based techniques can yield improved results when the model matches the true
conditions. The near-field propagation model, being based on an omnidirectional point

source in a homogeneous medium, does not perform well to localize the intra-thoracic

acoustic sources that were generated using an f; directional source. The transfer-function

based approach, being based on the true viscoelastic forward model, is able to provide
very accurate results. However, firrther studies need to be conducted to assess their
performance in cases of model-mismatch. The distributed source model approach
provides a scheme to improve the performance of the near-field propagation model and

may be used in absence of a transfer function model.

137

CHAPTER 6

SUMMARY AND CONCLUSIONS

6.1 Summary

A 2D finite 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 source-to-receiver
transfer functions corresponding to the various auscultation sites on the thorax are briefly
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-
field propagation model are studied and implemented. The localizing ability of these
algorithms to detect an intra-thoracic source are assessed and compared. In addition, four
methods based on a transfer-function approach are proposed to alleviate the issues of the
near-field propagation model based techniques. The results of the transfer fImction based
algorithms are compared to those obtained by the near-field 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-field

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 significant diagnostic information

o relation of the sensor signals to the extent and location of the damaged tissue or

the diseased state

0 type of multi-sensor 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 intra-thoracic 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), 974-987 (1997).

[2] A. Sarvazyan, “Audible-fiequency medical diagnostic methods: Past, present and
future,” J. Acoust. Soc. Am. 117, 2586 (2005).

[3] Y. Bahadirlar and H. O. Giilcfir, “Cardiac passive acoustic localization:
CARDIOPAL,” Elektrik, Turkish J. of Electrical Engineering & Computer Science 6(3),
243-259 (1998).

[4] N. L. Owsley and A. J. Hull, “Beamforrned nearfield imaging of a simulated coronary
artery containing a stenosis,” IEEE Trans. Medical Imaging 17(6), 900-909 (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, 1931-1946 (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), 428-436 (1963).

[8] J. J. Faber, “Damping of sound on the chest surface,” Circulation Res. 13, 352-358
(1963).

[9] L. G. Durand, J. Genest, and R. Guardo, “Modeling of the transfer function of the
heart-thorax acoustic system in dogs,” IEEE Trans. Biomed. Eng. BME-32(8), 592-601
(1985).

[10] J. Verburg, “Transmission of vibrations of the heart to the chestwall,” Adv.
Cardiovasc. Phys. 5(III), 84-103 (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, 248-249 (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, 925-934
(1989)

142

[13] G. R. Wodicka, A. Aguirre, P. D. DeFrain, and D. C. Shannon, “Phase delay of
pulmonary acoustic transmission fi'om trachea to chest wall,” IEEE Trans. Biomed. Eng.
39(10), 1053-1059 (1992).

[14] G. R. Wodicka, P. D. DeFrain, and S. S. Krarnan, “Bilateral asymmetry of
respiratory acoustic transmission,” Med. Biol. Eng. Comput. 32(5), 489-494 (1994).

[15] S. Lu, P. Doerschuk, and G. Wodicka, "Parametric phase-delay estimation of sound
transmitted through intact human lung," Med. Biol. Eng. Comput. 33(3), 293-298 (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, 667-676 (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, 1330-1335 (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), 25-31 (1996)

[20] Q. Grimal, A. Watzky, and S. Naili, “A one-dimensional model for the propagation
of transient pressure waves through the lung,” J. Biomechanics 35(8), 1081-1089 (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, 1109-1121 (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), 657-671 (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), 177-192 (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), 603-640 (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, 1029-1030 (1995).

[27] James V. Candy, “Model-based signal processing in acoustics: An overview,” J.
Acoust. Soc. Am. 101(5), 3155 (1997).

[28] Y. Bahadirlar and H. O. Gulcur, “Time-fiequency cardiac passive acoustic
localization,” Proceedings of the 23"! Annual International Conference of the IEEE
Engineering in Medicine and Biology Society, Vol. 2, 1850-1853 (2001).

[29] M. Kompis and G. R. ' Wodicka, “Spatial information in simultaneous
multimicrophone recordings of thoracic sounds,” J. Acoust. Soc. Am. 99(4), 2478-2500
(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, 1661-1664 (1998).

[31] M. Kompis, H. Pasterkamp, and G. R. Wodicka, “Acoustic imaging of the human
chest,” Chest 120, 1309-1321 (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, 117-122 (2005).

[34] A. M. McKee and R. A. Goubran, “Chest sound pick-up using a multisensor array,”
Proceedings of IEEE Conference on Sensors, 780-783 (2005).

[3 5] John Scmmlow and Ketaki Rahalkar, “Acoustic Detection of Coronary Artery
Disease,” Annual Review of Biomedical Engineering 9, 449-469 (2007).

[36] L. E. Kinsler and A. R. Frey, Fundamentals of acoustics (Wiley, New York, 3rd
edition, 1982).

[37] B. A. Auld, Acoustic fields 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, 707-714 (1951).

[39] A. Hosokawa, “Ultrasonic pulse waves in cancellous bone analyzed by finite-
difference time-domain methods,” Ultrasonics 44(1), e227-e23l (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/VI-I/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 soft-tissue conducted sound originating fi'om vocal
tract,” Applied Acoustics 70(3), 469-472 (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), 3665-3677 (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, 199-230 (2005).

[45] A. Hosokawa, “Simulation of ultrasound propagation through bovine cancellous
bone using elastic and Biot’s finite-difference time-domain methods,” J. Acoust. Soc.
Am. 118, 1782-1789 (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 low-frequency harmonic vibration,” Proceedings of the IEEE, Vol.
91(10), 1503-1519 (2003).

[48] X. Zhang, T. J. Royston, H. A. Mansy, and R. H. Sandler, “Radiation impedance of
a finite circular piston on a viscoelastic half-space with application to medical diagnosis,”
J. Acoust. Soc. Am. 109, 795-802 (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), 1346-1355 (1983).

[51] R. Natarajan, J. Williams, and G. Andersson, “Modeling changes in intervertebral
disc mechanics with degeneration,” J. Bone Joint Surgery Am. 88(2), 36-40 (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), 980-989
(2004).

[54] M. J ahed and S. J. Lai-Fook, “Stress wave velocity measured in intact pig lungs with
cross-spectral analysis,” J. Appl. Physiol. 76(2), 565-571 (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, 185-200 (1994).

[57] F. H. Drossaert and A. Giannopoulos, “A non-split complex fi'equency shifted PML
based on recursive integration for FDTD modeling of elastic waves,” Geophysics 72(2),
T9-Tl7 (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, 294-307 (2001).

[59] Q. H. Liu and J. Tao, “The perfectly matched layer for acoustic waves in absorptive
media,” J. Acoust. Soc. Am. 102(4), 2072-2082 (1997).

[60] J. Virieux, “P-SV wave propagation in heterogeneous media: velocity—stress finite-
difference method,” Geophysics 51, 889-901 (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), 474-481
(2002)

[63] G. Eason, J. Fulton, and I. N. Sneddon, “The generation of waves in an infinite
elastic solid by variable body forces,” Phil. Trans. Royal Soc. London (A) 248, 575-607
(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 finite-difference 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, 520-525 (2002).

[67] D. Lichenstein, G. Meziere, P. Biderrnan, and A. Gepner, “The 'lung point': an
ultrasound sign specific to pneumothorax,” Intensive Care Med. 29, 1434-1440 (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,
526-532 (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), 4-24 (1988).

[71] I. A. McCowan, D. C. Moore, and S. Sridharan, “Near-field Adaptive Beamformer
for Robust Speech Recognition,” Digital Signal Processing 12(1), 87-106 (2002).

[72] H. Krim and M. Viberg, “Two decades of array signal processing research: the
parametric approach,” IEEE Signal Processing Magazine 13(4), 67-94 (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), 1509-1522 (2002).

[74] Y. R. Zheng, R. A. Goubran, and M. El-Tanany, “Robust near-field adaptive
beamforrning with distance discrimination,” IEEE Trans. Speech and Audio Processing
12(5), 478-488 (2004).

[75] Y. Jin and B. Friedlander, “Detection of distributed sources using sensor arrays,”
IEEE Trans. Signal Processing 52(6), 1537-1548(2004).

147

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

IIIIII[III][IIIIIIIIIIIIIIIJIIIIIIIIIIII