.\ 73...}.

.2», $25? an...
..

 

 

 

.,
s n...
'- 1“}.

 

.v . .
.l f
.mmmfiunwflfifimfi
4r urémn... easifififlk
”unfit, “- Hfilfi
.JMdhfln. ,

 

 

 

 

Ilv
o.

b
. :uuwmuznmfi

 

 

 

 

 

‘ a , _ $35; _ ‘ :3;

 

 

 

mas-rs

:5; ,3 LIBRARY

" I Michigan State
University

 

 

This is to certify that the
thesis entitled

COMPARISON OF UWB SHORT-PULSE AND STEPPED-
FREQUENCY SYSTEMS FOR IMAGING THROUGH
BARRIERS

presented by

Benjamin Reid Crowgey

has been accepted towards fulfillment
of the requirements for the

- MS. degree in Electrical Engineering

 

WSW

Major Messor’s Signature‘

’3 'NOJCMImc " 2-009

 

Date

MSU is an Affirmative Action/Equal Opportunity Employer

—-—-—.—-—-
-——-—o
‘-—--
"o--—-
-.—-_.
—.—.--—-—-—. _

- ‘-_--.__
—-—---—-_
-—------
-~—-_- -

- ----_

PLACE IN RETURN Box to remove this checkout from yourrecord.
To AVOID FINES return on or before date due.

MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

SIDS K:/Prolecc&Pres/CIRC/DateDue.indd

COMPARISON OF UWB SHORT-PULSE AND STEPPED-FREQUEN CY
SYSTEMS FOR IMAGING THROUGH BARRIERS

By

Benjamin Reid Crowgey

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER SCIENCE
Electrical Engineering
2009

ABSTRACT

COMPARISON OF UWB SHORT-PULSE AND STEPPED-FREQUEN CY
SYSTEMS FOR IMAGING THROUGH BARRIERS

By

Benjamin Reid Crowgey

Ultra-wideband radar is useful for generating images of objects behind barriers due
to its penetrability and its capability to resolve small targets. Imaging systems may be
constructed to work directly in the time domain by generating baseband or monocycle
pulses and digitizing the received signal using a high-speed analog-to—digital converter,
or in the frequency domain by using a stepped-frequency narrow-band transmitter
and receiver. Each type of system has advantages and drawbacks. The time-domain
method provides ultra-wideband data with a single measurement, but it requires high-
speed A/ D converters to provide sufficient resolution, and adequate signal-to—noise
ratios may necessitate repetitive measurements. Stepped-frequency systems have high
dynamic range, but often require long dwell times and are subject to aliasing effects.

A canonical problem is established that allows the affects of various radar param-
eters on radar performance to be studied. For simplicity, a two-dimensional problem
is considered, consisting of a perfectly conducting strip located behind a lossy di-
electric slab of infinite extent illuminated by line sources. To assess the impact of
the parameters on system performance, images of the target are created using the
reflected field computed at several positions in front of the barrier and adjacent to
the sources. Specific parameters that are considered include sample rate, A/ D bit
length, pulse width, and SNR for a time-domain system, and sample rate, A/ D bit
length, bandwidth, and SN R for a frequency domain system. A time-domain labora-
tory system was constructed to investigate whether the image techniques used with

simulated data in the parameter study can be replicated in practice.

For Grandad

iii

ACKNOWLEDGMENTS

I would like to acknowledge my parents John and Martha Crowgey, who encour-
aged me to pursue my interests in engineering and helped me along the way. A special
thanks to my mother, who calmed me down when I needed it and was always there
to talk. I acknowledge my grandparents, Willard and Jenett Crowgey and Donald
and Jean Moore for their thoughts and support. Thank you to my sisters Sarah and
Julie Crowgey for their insight and encouragement.

I would like to thank Dr. Lee Kempel for taking a chance on me and showing
me the fun of Electromagnetics. If it wasn’t for the simple cell phone seminar six
years ago, who knows if I would even be here today. Thank you to Dr. Shanker
Balasubramaniam for your entertaining electromagnetic classes and support. Most
importantly, I would like to acknowledge my advisor, Dr. Edward Rothwell. Thank
you for all your help and support on this research. I know I wouldn’t have made it
through without your guidance and light humor. Even if you had some interesting
comments on my writing techniques, I know you always had my best interests at
heart.

I acknowledge Dr. Eric Mokole, head of the Surveillance Technology Branch,
Radar Division, at the Naval Research Laboratory in Washington, DC. for initiating
the working relationship and providing the funding for this research effort.

A special thanks to Olukorede A Akinlabi-Oladimeji for spending extra time doing
some very tedious measurements with me. Also, thank you to Collin Meierbachtol
for your random help throughout my masters career. You have both become great
friends.

I would like to acknowledge the Michigan State University Electromagnetics Re-
search Group for providing assistance on all of the projects I have worked on at MSU.

It has been an honor to work with so many talented people, many of whom will be

iv

lifetime friends.

Most of all I would like to thank my girlfriend, Stacey Richardson, for all her help
and support during these final months. You are a great friend and a terrific person. I
am lucky to have found someone like you. You spent hours proofreading this thesis,

and I know Dr. Rothwell and I both thank you.

Table of Contents

List of Tables .................. . ................ viii

List of Figures .................................. xv

1 Introduction ................................. 1

1.1 Background ................................ 2

1.2 Research Overview ............................ 2

2 Canonical Problem Theory ........................ 4
2.1 Electric Fields impressed by a current source in the Presents of a Di-

electric Slab ................................ 5

2.1.1 Current Source .......................... 6

2.1.2 Partial Differential Wave Equation Derivation ......... 6

2.1.3 Ordinary Differential Equation Derivation ........... 11

2.1.4 Relationship Between Vector Potentials and Fields ...... 14

2.1.5 Solutions to ODE ......................... 15

2.1.6 Electric Field Derivation ..................... 25

2.2 Scattered Field from Conducting Strip ................. 26

2.2.1 Integral Equation Derivation ................... 27

2.2.2 Amn Evaluation ......................... 31

2.2.2.1 Umn Evaluation .................... 33

2.2.2.2 an Evaluation .................... 37

2.2.3 Matrix Equation ......................... 40

3 Numerical Results ............................. 41

3.1 Computation of Incident and Scattered Fields ............. 42

3.2 Imaging Analysis ............................. 44

3.2.1 Frequency to Time Domain Transformation .......... 44

3.2.2 2D Image Construction ...................... 48

3.2.3 Image Quantification ....................... 51

4 Parameter Analysis ............................. 54

4.1 Time Domain Analysis .......................... 54

4.1.1 Pulse Width ............................ 54

4.1.2 Sampling Rate .......................... 65

4.1.3 Signal-to-Noise Ratio ....................... 76

vi

4.1.4 Digitization ............................ 87

4.2 Frequency Domain Analysis ....................... 98
4.2.1 Bandwidth ............................. 99

4.2.2 Sampling Rate .......................... 109

4.2.3 Signal-to—Noise Ratio ....................... 116

4.2.4 Digitization ............................ 127

4.3 Two Conducting Strips .......................... 139
4.3.1 Pulse Width ............................ 141

4.3.2 Bandwidth ............................. 151

5 Experiment .................................. 163
5.1 Experiment Set up ............................ 164
5.2 Experiment Calibration .......................... 169
5.3 Experiment Results ............................ 173

6 Conclusion .................................. 175
6.1 Topics for further study ......................... 176
Appendix A: Fortran Code for Electric Scattered Field ........ 178
Appendix B: Fortran Code for Electric Incident Field ......... 201
Appendix C: Matlab Imaging Code .................... 210
Bibliography ................................... 225

vii

List of Tables

6.1 Minimum system requirements for a time domain radar system. . . . 176

6.2 Minimum system requirements for a frequency domain radar system. 176

viii

List of Figures

2.1

2.2

2.3

3.1

3.2

3.3

3.4

3.5

3.6

3.7

3.8

3.9

4.1

4.2

4.3

4.4

4.5

Canonical problem geometry. ...................... 4
Partitioning of canonical problem geometry. .............. 5
Conducting strip parameters. ...................... 27
Position of source/ observation pairs. .................. 41

Magnitude of the frequency domain scattered field at a position cen-
tered above the strip ............................ 45

Magnitude of the frequency domain scattered field weighted with a GMC. 46

Time-domain scattered field at a point centered above the conductive

strip ..................................... 47
Time-domain incident field at a point centered above the conducting

strip ..................................... 48
Image of lossy dielectric slab and conducting strip.. .......... 50

Proportionally correct mage of lossy dielectric slab and conducting strip. 50
Area in which sharpness of conducting strip is determined. ...... 52

Circle of radius equal to the image radius centered around the con-
ducting strip with scale x1011 ....................... 53

Effect of pulse width on image radius ................... 55
Conducting strip image constructed using a pulse width of T = 0.06 ns. 56
Conducting strip image constructed using a pulse width of T = 0.10 ns. 57
Conducting strip image constructed using a pulse width of T = 0.16 ns. 58

Conducting strip image constructed using a pulse width of T = 0.26 ns. 59

4.6

4.7

4.8

4.9

4.10

4.11

4.12

4.13

4.14

4.15

4.16

4.17

4.18

4.19

4.20

4.21

4.22

4.23

4.24

Conducting strip image constructed using a pulse width of T = 0.50 ns. 60
Conducting strip image constructed using a pulse width of T = 1.00 ns. 61
Conducting strip image constructed using a pulse width of T = 1.50 ns. 62

Conducting strip image constructed using a pulse width of T = 2.00 ns. 63

Conducting strip image constructed using a pulse width of T = 2.50 ns. 64
Conducting strip image constructed using a pulse width of T = 3.00 ns. 65
Effect of sampling period on image radius. ............... 66
Conducting strip image constructed using a sampling period of At =
0.05 ns. .................................. 67
Conducting strip image constructed using a sampling period of At =
0.10 ns. .................................. 68
Conducting strip image constructed using a sampling period of At =
0.20 ns. .................................. 69
Conducting strip image constructed using a sampling period of At =
0.30 ns. .................................. 70
Conducting strip image constructed using a sampling period of At =
0.40 ns. .................................. 71
Conducting strip image constructed using a sampling period of At =
0.50 ns. .................................. 72
Conducting strip image constructed using a sampling period of At =
0.60 ns. .................................. 73
Conducting strip image constructed using a sampling period of At =
0.70 ns. .................................. 74
Conducting strip image constructed using a sampling period of At :
0.80 ns. .................................. 75
Conducting strip image constructed using a sampling period of At =
1.60 ns. .................................. 76
Effect of noise on image radius. ..................... 77
Conducting strip image constructed using a SN R = 20 dB. ...... 78

4.25 Conducting strip image constructed using a SN R = 18 dB. ...... 79

4.26 Conducting strip image constructed using a SNR = 16 dB. ...... 80
4.27 Conducting strip image constructed using a SNR = 14 dB. ...... 81
4.28 Conducting strip image constructed using a SNR = 12 dB. ...... 82
4.29 Conducting strip image constructed using a SNR = 10 dB. ...... 83
4.30 Conducting strip image constructed using a SNR = 8 dB ........ 84
4.31 Conducting strip image constructed using a SNR = 6 dB ........ 85
4.32 Conducting strip image constructed using a SNR = 4 dB ........ 86
4.33 Conducting strip image constructed using a SNR = 2 dB ........ 87
4.34 Effect of digitization on image radius ................... 88
4.35 Conducting strip image constructed using the number of bits 2 12. . 89
4.36 Conducting strip image constructed using the number of bits = 11. . 90
4.37 Conducting strip image constructed using the number of bits = 10. . 91
4.38 Conducting strip image constructed using the number of bits = 9. . . 92
4.39 Conducting strip image constructed using the number of bits I: 8. . . 93
4.40 Conducting strip image constructed using the number of bits = 7. . . 94
4.41 Conducting strip image constructed using the number of bits = 6. . . 95
4.42 Conducting strip image constructed using the number of bits = 5. . . 96
4.43 Conducting strip image constructed using the number of bits 2 4. . . 97
4.44 Conducting strip image constructed using the number of bits = 3. . . 98
4.45 Effect of digitization on image radius ................... 99

4.46 Conducting strip image constructed using a bandwidth = 110.38%. . 100

4.47 Conducting strip image constructed using a bandwidth 2 36.75%. . . 101
4.48 Conducting strip image constructed using a bandwidth = 22.13%. . . 102
4.49 Conducting strip image constructed using a bandwidth 2 15.76%. . . 103

xi

4.50 Conducting strip image constructed using a bandwidth = 12.26%. . . 104

4.51 Conducting strip image constructed using a bandwidth = 10.03%. . . 105
4.52 Conducting strip image constructed using a bandwidth 2 8.49%. . . . 106
4.53 Conducting strip image constructed using a bandwidth = 7.35%. . . . 107
4.54 Conducting strip image constructed using a bandwidth = 6.49%. . . . 108
4.55 Conducting strip image constructed using a bandwidth = 5.81%. . . . 109
4.56 Effect of sampling rate on image radius. ................ 110
4.57 Conducting strip image constructed using a sampling rate of A f = 20

MHz ..................................... 111
4.58 Conducting strip image constructed using a sampling rate of A f = 40

MHz ..................................... 112
4.59 Conducting strip image constructed using a sampling rate of A f = 60

MHz ..................................... 113
4.60 Conducting strip image constructed using a sampling rate of A f = 80

MHz ..................................... 114
4.61 Conducting strip image constructed using a sampling rate of A f = 100

MHz ..................................... 115
4.62 Conducting strip image constructed using a sampling rate of A f = 120

MHz ..................................... 116
4.63 Effect of the SN R on image radius. ................... 117
4.64 Conducting strip image constructed using a SN R = 44 dB. ...... 118
4.65 Conducting strip image constructed using a SNR = 42 dB. ...... 119
4.66 Conducting strip image constructed using a SNR = 40 dB. ...... 120
4.67 Conducting strip image constructed using a SNR = 38 dB. ...... 121
4.68 Conducting strip image constructed using a SNR = 36 dB. ...... 122
4.69 Conducting strip image constructed using a SNR = 34 dB. ...... 123
4.70 Conducting strip image constructed using a SNR = 32 dB. ...... 124

xii

4.71 Conducting strip image constructed using a SN R = 30 dB. ...... 125

4.72 Conducting strip image constructed using a SNR = 28 dB. ...... 126
4.73 Conducting strip image constructed using a SNR = 26 dB. ...... 127
4.74 Effect of digitization on image radius ................... 128

4.75 Conducting strip image constructed using the number of bits 2 12. . 129
4.76 Conducting strip image constructed using the number of bits 2 11. . 130

4.77 Conducting strip image constructed using the number of bits = 10. . 131

4.78 Conducting strip image constructed using the number of bits = 9. . . 132
4.79 Conducting strip image constructed using the number of bits = 8. . . 133
4.80 Conducting strip image constructed using the number of bits = 7. . . 134
4.81 Conducting strip image constructed using the number of bits 2 6. . . 135
4.82 Conducting strip image constructed using the number of bits 2 5. . . 136
4.83 Conducting strip image constructed using the number of bits = 4. . . 137
4.84 Conducting strip image constructed using the number of bits 2 3. . . 138
4.85 Conducting strip image constructed using the number of Bits = 2. . . 139
4.86 Parameters of two conducting strip study ................ 140
4.87 Effect of pulse width on image radius of two conducting strips ..... 141

4.88 Two conducting strip image constructed using a pulse width of T =
0.05 ns ................................... 142

4.89 Two conducting strip image constructed using a pulse width of T =
0.10 ns ................................... 143

4.90 Two conducting strip image constructed using a pulse width of T =
0.15 ns ................................... 144

4.91 Two conducting strip image constructed using a pulse width of T =
0.25 ns ................................... 145

4.92 Two conducting strip image constructed using a pulse width of T =
0.50 ns ................................... 146

xiii

4.93 Two conducting strip image constructed using a pulse width of T =
1.00 ns ...................................

4.94 Two conducting strip image constructed using a pulse width of T =
1.50 ns ...................................

4.95 Two conducting strip image constructed using a pulse width of T =
2.00 ns ...................................

4.96 Two conducting strip image constructed using a pulse width of T = 2.5
ns .....................................

4.97 Two conducting strip image constructed using a pulse width of T =
3.00 ns ...................................

4.98 Effect of bandwidth on image radius of two conducting strips.

4.99 Two conducting strip image constructed using a bandwidth of 110.38%

4.100Two conducting strip image constructed using a bandwidth of 36.75%
4.101Two conducting strip image constructed using a bandwidth of 22.13%
4.102Two conducting strip image constructed using a bandwidth of 15.76%
4.103Two conducting strip image constructed using a bandwidth of 12.26%
4.104Two conducting strip image constructed using a bandwidth of 10.03%
4.105Two conducting strip image constructed using a bandwidth of 8.49%
4.106Two conducting strip image constructed using a bandwidth of 7.38%
4.107Two conducting strip image constructed using a bandwidth of 6.49%

4.108Two conducting strip image constructed using a bandwidth of 5.81%

5.1 Geometry of arch-range scattering system at MSU. ..........
5.2 Time domain scattering measurement configuration ...........
5.3 Pulse generated by the PSPL 5208-DC pulse generating network.

5.4 Spectrum of the generated pulse. ....................

5.5 Location of transmitting and receiving antennas .............

xiv

147

148

149

151

152

153

154

155

156

157

158

159

160

161

162

164

165

166

168

5.6 Experimental setup. ........................... 169

5.7 Block diagram model of measurement system. ............. 170
5.8 Raw measured waveform for a conducting plate. ............ 172
5.9 System transfer function .......................... 173

5.10 Image of barrier and strip created using time-domain response at re-
sponses measured at 21 different transmitter/ receiver positions. . . . 174

XV

Images in this thesis are presented in color

 

xvi

Chapter 1

Introduction

A comparison between ultra-wide band (UWB) short-pulse and stepped—frequency
systems for imaging through a barrier is undertaken in this thesis. The penetration
ability of UWB radar, and its capabilities for resolving small targets, makes this tech-
nology quite appealing for producing images of objects behind barriers. The concept
of a UWB radar can be used to implement both time and frequency domain systems.
A time-domain—sytem may be implemented by generating baseband or monocycle
pulses and then digitizing the received signal with the use of a high—speed analog-to-
digital (A/ D) converter. In a frequency domain system, a stepped-frequency narrow-
band transmitter and receiver may be used. With each type of system there are
advantages and disadvantages. For the time-domain method the UWB signal pro-
vides data with a single measurement, but requires high-speed A/ D converters to
ensure sufficient resolution of the image. In addition, inadequate signal-to—noise ratio
(SNR) may produce a need for repetitive measurements and averaging, increasing the
acquisition time. In the frequency domain, long dwell time and aliasing effects are

concerns, while high dynamic range is a significant advantage.

1. 1 Background

A significant amount of research has been done to analyze scattering from objects
located behind barriers [1]. In addition, various processing methods have been de-
veloped to obtain images of the scatterer [2]. This thesis focuses on the comparison
between short-pulse and stepped-frequency systems for imaging through a barrier [3]
[4] [5]. The UWB simulations or measurements are performed in a synthetic aper—
ture radar (SAR) configuration [6] [7] [8] and then images are created using a simple

scattering—center technique.

1 .2 Research Overview

A two-dimensional canonical problem is established to study the affects of various
parameters on the quality of the image produced by time-domain and stepped fre-
quency UWB radar systems. This canonical problem consists of a line source located
above a lossy dielectric slab representing a barrier. The interrogating cylindrical
waves produced by the line source impinge on the barrier and interact with targets
placed behind the barrier. The scattered field is computed in the frequency domain
at specific discrete frequencies within a chosen band. This simulates the data that
would be acquired by a UWB stepped-frequency system. An equivalent time-domain
system is studied by using an inverse Fourier transform to convert the data into the
time domain.

Images are constructed using a simple scattering—center technique to provide a
means of assessing radar performance. Factors such as pulse width, sampling rate,
number of bits, SN R, and time jitter are examined in the time domain to determine
their affect on the quality of the image as measured in terms of the dispersion of
localized scattering centers. Bandwidth, sampling rate, SNR, and number of bits are

examined in the frequency domain.

In addition, a laboratory time—domain radar system was constructed using instru—
mentation to allow a validation of the simulation results. The laboratory system is
not a dedicated radar, but rather uses standard instrumentation. The data was ac-
quired using the Michigan State University (MSU) reflectivity arch range consisting
of two horn antennas. A barrier was constructed to replicate the simulated canonical

problem.

Chapter 2

Canonical Problem Theory

To investigate the different properties of ultra-wide band short pulse and step-frequency
imaging radar systems, a simple canonical problem is considered. Using simulations
instead of measured data allows the affects of important system parameters to be
studied over a much broader range of values than is allowed through laboratory ex—
perimentation. The canonical problem, shown in Figure 2.1, consists of an infinite
electric line source aligned along the x-direction located at a height 2 = h above a

lossy dielectric slab, which represents a barrier. This dielectric slab has a thickness t,

[Line Source Lossy Planar z=h
Barrier

£0 1 Q g z-
z _
/ 8,0’ -0

8 _
0 y Z=Zs-

 

 

 

 

 

 

Figure 2.1: Canonical problem geometry.

and is infinite in the x- and y-directions. The slab has relative permittivity er, con-

ductivity or, and free—space permeability no. A target consisting of an infinitesimally
thin perfectly-conducting (PEC) strip is located a distance z = 23 behind the slab.
To solve for the electric fields for this canonical problem, the field incident on the
conducting strip in the presence of the slab must first be determined. The boundary
condition of zero total tangential electric field on the surface of the strip may then
be employed to derive an integral equation for the surface current induced on the
strip. This current is then used to compute the scattered field, and then the total
field, which is what is measured in practice, is found by summing the incident and

scattered fields.

2.1 Electric Fields impressed by a current source
in the Presents of a Dielectric Slab

The fields incident on the conducting strip, in the presence of a dielectric slab, are
derived by first considering Figure 2.2, where the problem is separated into 4 different

regions. Regions 2 and 3 are chosen to be two separate regions because this results

 

 

Region 3

——————————— ©—-—--—-——-— -———— Z=h
Region 2 8
Egan: Z//////////?Y/y//////////// e ,0 :3,
Region 0 80

Figure 2.2: Partitioning of canonical problem geometry.

in source free regions with the line source located on an interface.

2.1. 1 Current Source

The infinite line source is positioned at z = h and y = 0; therefore the volume current

density can be written as

—o

J(y, z) = "16(z — h)6(y). (2.1)

The surface current density can be determined through the use of boundary conditions

in this 2D canonical problem [9]. This current is given by

h+A

K(y)=Ali:I+10/J(y,z)dz. (2.2)

Substituting the volume current density gives

h+A

my) = Aline i16(z—h)6(y)dz. (2.3)

The surface current density for the infinite line source is thus

Re) = My). (24)

2.1.2 Partial Differential Wave Equation Derivation

The fields are represented in terms of potentials to facilitate their computation. The
magnetic flux density vector, B, is expressed as the curl of the magnetic vector po-
tential, A:

ézvXA am

Substituting (2.5) into the differential form of Faraday’s law,

- as
V x E — _—df’ (2.6)
results in
- at

The sum of the two vector quantities can thus be expressed as the gradient of the

electric potential, (15:

- at
E + E — —V¢. (2.8)
Thus,
- at
E = —V ' — —. 2.
<0 8t ( 9)
The differential form of Gauss’s law is
V - D = p (2.10)

or, from the constitutive relations assuming a homogeneous medium,

VB: ’3. (2.11)

6

Substituting (2.9) into (2.11) results in

821‘
01‘
_ 2 _ 2 . -' _ e
V at (V A) _ 6, (2.13)

which is a differential equation for (13 and A.

Next, examine the differential form of Ampre’s law,

.. .. 313
V H = —.
x J + (9 t
The constitutive relation,
§=pfi,

along with (2.5) is substituted into (2.14) to give

—o

VXVXAzpj+p%lt2.

an)

an)

(2.16)

In a conducting medium, as seen in this canonical problem, the source field must

be separated into two parts. The first is a causative impressed term j;- which is

independent of the produced fields. Also necessary is a secondary term is which

exists only as an effect of the sourced fields. For an isotropic conducting medium the

effect of the sourced field is given by Ohm’s law,

3:03

Writing the total current as

(2.16) becomes,

VXVXAzpaE+p%—It)+j%.

Then, using the constitutive relation

(2.17)

(2.18)

(2.19)

am)

(2.19) is rewritten as
V XVX/I= (pa+pe—)E+fi.
Substituting (2.9) into (2.21) yields
a a 621' -
V XV XA— (INTI—HEB?) (—V¢—-(737) +Jz

or
.. a 621' 02/?
V x V X A = —,uaV¢ — pean) — [Ia—a- — [16 at?

 

Recalling the vector identity

—o —o

VxVxA=V(V-A) v2.4,

(2.23) can be rewritten as

at 824 -

V2A=V(V-A+pa¢+pe%) +p0—+u6—8—t§—JZ~.

(9t at

To simplify (2.25), the Lorentz condition for potentials

18¢_

.A'=———
V V2 at

limb

is used, where z/ is given by

 

This gives

V - A = fine? —/10(,73.

+.7,-.

(2.21)

(2.22)

(2.23)

(2.24)

(2.25)

(2.26)

(2.27)

(2.28)

Inserting (2.28) into (2.25), some terms on the right side of (2.25) vanish, thus giving

62/? -

V2A= pa— 6A e— 3t2 —J,. (2.29)

at +“

This is a time domain wave equation for the vector potential A By substituting

(2.28) into (2.13) results in

2
—V2$—pc2—‘2—,m g—f— f,

(3t2 (2.30)

which is a wave equation for the scalar potential.
The analysis performed in the remainder of this thesis will be undertaken in the
frequency domain. Thus, a frequency domain wave equation is required. Transform-

ing (2.29) into the frequency domain through the Fourier transform identity

gt (a) jw, (2.31)
gives
V2.4 = no (jw) 21+ ,ue (2 (.12) xi— J, (2.32)
or
372/1. 2 (jump — are?) A— .12. (2.33)

All quantities are new frequency domain variables. Define the wavenumber as
k = (32,166, (2.34)
where the complex permittivity is

e =e—j—. (2.35)

10

Then (2.33) reduces to
v24 + T34 = —.f,. (2.36)

This is a partial differntial equation for the vector potential, which is solved using in
subsequent sections using spatial Fourier transform techniques.
As shown in Figure 2.2, the geometry is partitioned into 4 regions. Region 1,

inside the dielectric slab, has the wavenumber, k, defined through
k2 = (422/1060. (2.37)
Since this is a source free region, (2.36) becomes
v24 + 162/T = 0. (2.38)
In regions 0, 2, and 3 the wavenumber is given through
168 = (32,1060, (2.39)

because these regions are free space. Also, these regions are source free. Therefore

the partial differential wave equation (2.36) becomes

v24 + 163.4 = 0. (2.40)

2.1.3 Ordinary Differential Equation Derivation

An ordinary differential equation (ODE) is a relation that contains functions of only
one independent variable. While a partial differential equation (PDE) is a relation
which contains partial derivatives of several variables. Thus, ODEs are computation-

ally easier to solve and are homogeneous, for this canonical problem, since the line

11

source lies on the interface between region 3 and region 2. To reduce the PDEs to
ODEs, it is noticed the source and fields are invariant in the sis-direction. Therefore a
Fourier transform of the y variable is performed to reduce the PDE to an ODE. The

definition of the Fourier transform pair is

00

A(ky,z) = /A(y,z)e_jkyydy (2.41)
—00
1 oo

A(y,z) = 2—7; / A(ky,z)e]kyydky. (2.42)

Inserting (2.42) into the partial differential wave equation for region 1, (2.38), results
in

00
(V2+k2) )2i /( 4(a) )eJ flaky: 0 (2.43)
—00

Since there is only a a: component of the current seen in Figure 2.2, the vector magnetic

potential A only has a a: component:
4(a), z) = 2.433067), 2). (2.44)

Therefore (2.43) is rewritten as

00
(V2+k2) 71” / A$(ky,z)ejk3/ydky :0. (2.45)
”00
Ol' 00
a2 a2 a? 2 1 - 4,. .
__ __ —— k: e} yydk. = . 2.4
(3$2+8y2+02——2 +1“ 26 f A“ 3”") y 0 ( 6)

—OO

12

Since the integral is only dependent on y and z, the partial derivative with respect

to x is zero. Performing the remaining derivatives results in

00 OO
1 ~ ° ~ 1 ~ - .
[€2- / AQIUCy, Z)€Jkyydky — 2—7? / A1;(ky, Z)k§€]kyydky 'I'

—oo —oo

00 ..
1 AHA/971,2) jkf,
— —‘ ‘Jydk = .
271’ 832 6 y 0

Factoring out the common variables from this expression, (2.47) reduces to

27r
—oo

00
1 2 (92 2 ~ "A:

The Fourier integral theorem then implies

02 ~
(52:2— +132) Ax(k‘y, Z) = 0,

with
p = ilfkg — 163

in regions 0, 2, and 3. For region 1, (2.48) requires

62 ~
(5;? + 02> “(2’2 z 0’

where

_ ,/ 2_ 2
q—i k kg.

(2.47)

(2.48)

(2.49)

(2.50)

(2.51)

(2.52)

Note that the proper choice of signs on the radicals depends on the physics of the

problem and is considered later.

13

2.1.4 Relationship Between Vector Potentials and Fields

The electric and magnetic fields can be found from the vector potentials once the

ODES (2.49) and (2.51) are solved. To find this relationship first recall the divergence

of A:
~ 8A$ aAy 8A3
V - A = —— — —. 2.53
8:1: + By + 32 ( )
Since there is only an :1: component of A, the divergence reduces to
-' 8A$(y~ 2)
V - A = ———’—. 2. 4
0x ( 5 )
But A3; is only dependent on y and 2, so
V 4:0 asm

Writing the Lorentz condition in the frequency domain and using (2.55) thus gives

3:9 (2%)

The electric field is then, from (2.9),

E = —ij. (2.57)

Combining (2.5) and the constitutive relation (2.15), gives the magnetic field as

1 _.

 

 

Hz—VXA, 9%)
,u
or from the definition of the curl,
_. 1 19/151; .8A3;
H = — — . 2.59
M (y 32 Z By ) ( I

14

The relations between the fields and the vector potentials are thus summarized as

1 8A3;
H == —— 2.61
y #0 az ( )
Hz 2 __L2_A_£ (2.62)
#0 33/

2.1.5 Solutions to ODE

The ordinary differential equations ( 2.49) and (2.51) are each the harmonic differential
equation. Thus both have solutions that are complex exponential or trigonometric

functions. In region 3, z > h, there is only an upward propagating wave. Therefore

~

A1;(k1,z) = Ole—4P2, (2.63)

where 01 is a constant to be determined using appropriate boundary conditions. In
region 2, t 5 z < it, there are both upwards and downwards propagating waves,

requiring

~

Away, 2) = CQijZ + Cge—jpz. (2.64)
In region 1, 0 _<_ z < t, the potential can be written in terms of standing waves:

~

A5,;(ky, z) = C4 sin qz + C5 cos qz. (2.65)

Region 0 has only a downwards propagating wave. Therefore

~

Amy, z) = 06697”. (2.66)

It is important that the appropriate sign on p and q from equations (2.50) and (2.52)

be determined through the use of physical reasoning. In region 3, when k; 3 k2, the

15

wave must propagate upwards. Therefore, the sign on expressions in the exponent of

(2.63) must be negative. This requires

9 = +(/k(2 — k2, 6% _<_ kg. (2.67)

When 16% > 168, the wave must upward evanescent, again requiring the expression in

the exponent to be negative. Thus,

,0 = _j,/k.§ — 162, 1:5 > 1.3.9.68)

This is equivalent to requiring Re{p} Z 0 and Im{p} S 0. Through this same type
of physical reasoning, it can be determined that Re{q} Z 0 and Im{q} S 0.

The boundary conditions of the field vectors at the interfaces of different media are
necessary to solve electromagnetic canonical problems involving contiguous regions
of different parameters. These boundary conditions are found by applying Maxwell’s
equations to a small region at the interface of two different media. From these deriva-
tions it is found that the tangential component of the electric field is continuous across

an interface. In other words
a x (131 — E2) = 0, (2.69)

where ft is the normal vector pointing into region 1 from region 2 , t designates the
tangential component of E field and the numbers designate the region on each side of
the interface. The normal component of D is discontinuous across an interface where
a surface charge exists, with the amount of discontinuity being equal to the surface

charge density:

In contrast, the normal component of B is continuous across an interface and the

16

tangential component of H is discontinuous across an interface where a free surface

current exists. Thus,

ft- (81 -- 82) = 0 (2.71)

6 x (H1 — H2) = J3 (2.72)

With this set of boundary conditions, the constants in the solutions to the ODEs
can be found. First, the respective E fields are found from (2.60) and the boundary

conditions between regions 0 and 1 are employed. Since (2.69) holds true at z = 0,

[C4 sin zq + C5 cos zq]z=0 = [CGejpz] (2.73)

2:0,

and thus,
C5 = C6. (2.74)

The tangential H fields are examined next. Using (2.72) and (2.61) to derive H from

the potentials, and noting that there is no current source on the surface 2 = 0,

1 6 1 8 '4
—— C sinz +C cosz ] = [—— C 631)“ ] . 2.75
#0 5z( 4 q 5 (1) 2:0 #0 32( 6 ) 2:0 ( )

Taking the derivative with respect to z and multiplying through by #0 gives

 

 

(C4q cos zq — C5q sin zq) = C6jpejpz (2.76)
z=0 z=0
Inserting z = 0 thus produces
C4 = g 6' (2.77)

On the boundary between regions 1 and 2 E is continuous, resulting in

[_ij2ejpz — ij3e_ij] z—t = [—jw (C4 sin zq + C5 cos zq)]z=t. (2.78)

17

Substituting for C4 and C5 from (2.77) and (2.74), (2.78) becomes
[—ij2ejpz — ij3e—jPZ] z—t = [C6 (% sin zq — jw cos 2(1)] . (2.79)

Thus,
02,3th + C3e_jpt = C6 (353 sin zt + cos 2t) . (2.80)

Next, the boundary conditions for the H fields are employed at z = t. Using (2.72)

gives

[—1--3 (C26jpz + C3e_jpz)] = [95—9— (23 sin zq + cos zq)] . (2.81)
#0 az 2:7 #0 52 q z=t

Taking the derivatives with respect to z and substituting z = t results in
ijerpt — ij3e-jpt = C6 (jpcos zt — qsin qt). (2.82)
Continuity of E across the boundary between regions 2 and 3 requires

[—ij26jpz —jwe—jpz] = [—ijle_jpz] h (2.83)

2:

z=h

according to (2.69). Substituting 2 = h then gives
Cgejph + e_jph = Cle—jph. (2.84)

On the boundary between regions 2 and 3, the tangential H field is discontinuous by

the surface current, resulting in

2 x (11ng — Hgyg) : K352. (2.85)

18

Where If}; is the transform domain surface current density given, using (2.4), as

00

K~x(ky, z) = / Id(y)e_jkyydky. (2.86)
--00
Integrating gives
[6213(ky, z) : I. (2.87)
With this, (2.85) becomes
(9 1 _ - '3 1 ,- _ -
___c,. m _ __ (-8298 _ 63.. m) = 1. (2.88)
02 #0 5z #0 z=h
Differentiating and substituting z = h yields
—ijle_jph — ijerph +j12C3e_jph = —1,u0. (2.89)
In summary, there are 6 equations and 6 unknowns,
C5 = CG (2.90)
C4 = gCG (2.91)

Cgejpt + C3e_jpt = C6 (% sin zt + cos at) (2.92)

.7190 ijt — jPC 63—th = C jpcos zt — qsin qt (2.93
2 3 6

CQB-jph + e’jph : ole-3'1)" (2.94)

The 6 constants in (2.90)-(2.95) can be found using tedious standard linear algebra

19

techniques. First, (2.92) is multiplied by jp and added to (2.93), giving

. 2
2C2jpejpt = C6 (JD? sin qt + 2jpcos qt — qsin qt) , (2.96)
or
C6 192 . . .
C = —.— ——s1n t+2 cosqt — smqt . 2.97
2 23.106”, ( q q 37) q ) ( )

Then the multiplication of (2.92) with jp is subtracted from (2.93), resulting in

. _ 2
2C3jpe—th = C6 (Pq— sin qt + qsin qt) , (2.98)
or
C6 132 . .
C =———.—— —sm t+ smt . 2.99
3 ije_,p,(q q q q) ( )

Next, substitution of C2 and C3 from (2.97) and (2.99) into (2.94) gives

C -. _ 2
_.6_e]1)(h t) _I)_ sin qt + 2jp COS qt — qsin qt 'I'
231) (I

C _. _ 2 _,-

Multiplying through by 2 jp and factoring out C6, gives the equation

C6

 

. 2
er(h_t) (_p_ sin qt + 2jpcos qt — qsin qt) +
q

. 2 .
e—Jp(h—t)(p—- sin qt + qsin qt) = 2C1jpe-Jph. (2.101)
(I

 

To continue the procedure, C2 and C3 from (2.97) and (2.99), are substituted into

20

 

(2.95), yielding

C . __ 2
32831701 t) (_%— sin qt + 2jpcos qt — qsin qt)

. 2 .
_%e—JP(h—tl (Eq— sin qt + qsin qt) = 1,120 — Cljpejph.

Rearranging gives

06

 

. 2
(swat—'5) (—% sin qt + 2jpcos qt — qsin qt)

. 2 .
—e—Jp(h—t) (2;— sin qt + qsin qt)] = 21710 — 2C1jpejph.

To simplify notation, (2.101) is rewritten as

CG” [X1 01’

where

. 2
X = [e]p(h—t) (——8— sin qt + 2jpcos qt —— qsin qt)
9

. 2
677120143) (Eq— sin qt + q sin qt)
Similarly, (2.103) is written as

06 [Y] : —21p.0 + 2C1jpe—jph,

21

].

(2.102)

(2.103)

(2.104)

(2.105)

(2.106)

where

Y:

 

._ 2
— ejp(h_t) (—-PL sin qt + 2jpcos qt — qsin qt)
‘1

+
. h p2
e—Jp( —t) —q- sin qt + q sin qt . (2.107)

Substituting (2.104) into (2.106) yields

2;, —jph _ -.
[Y]-%k—]—C1 = —2[p.0 + 2C1jpe Jph, (2.108)
or ‘
2- —]])h _ .
([Y]J—pffi— — 2jpe Jph) C1 = —2[11.0. (2.109)

Factoring out similar constants results in

C12jpewjphCKI—[i71-E—I) = —2[1¢0, (2.110)
yielding
_ 1136' 11 [XI
01‘ 2'19 31) (in—[X1] (2'1“)

As seen in (2.111), [Y] — [X] must be computed. This results in the equation

. 2 2
Y — X = er(h_t) (1); sin qt + 3— sin qt — 2jpcos qt — 2jpcos qt + qsin qt
q

. _ 2 2
+q Sin qt) + e—]p(h—t) (_p_ sin qt + 1)— sin qt + qsin qt — qsin qt) (2.112)
q q

or, by subtracting out like terms,

' I t —p2
Y — X = —2e]p( If I — sin qt + 2jpcos qt — qsin qt . (2.113)
q

22

Therefore

IX] jp(h—t) P2 - - .
————= 6 -—sm t+2-cost— smt +
. h p2
e—]p( _t) —sinqt+qsinqt x
q
. _ 2 -1
[—2eJP(h‘t) (1 sin qt + 2jpcos qt — qsin 74)] . (2.114)
‘1

Simplifying (2.114) yields

2

 

 

[X] 1 e‘2jPIh—t) :42’ Si“ qt + 251“ qt (2 115)
—. = —— — - 2 . .
[Y X] 2 2 7% sin qt + 2J1? COS (It — q sin qt

Inserting (2.115), into (2.111) produces C1:

_22 , _
Cl = fliejph + lu—OeTQth-tlejl’h 2 q s1nqt+qs1nqt
ZJP 2]}? :qE' sin qt + 2jp cos qt — qsin qt

(2.116)
The next constant to be solved for is C6. Substituting C1 from (2.110) into (2.104),

a simplified equation for C6 results:

 

—2I;LO
= ——————. 2.117
06 [Y] — IX] ( )
Inserting (2.113) gives C6:
e-J'IIUI-t)
06:1/10— 2 . (2.118)

—qL sin qt + 2jp cos qt — qsin qt

23

 

Substituting C6 into (2.97) then gives C2:

C =flflle—JPUl—tle—th ‘1

2
2' __ 2
Jp —qB—sin qt + 2jpcos qt — qsin qt

This is simplified to give
[#0 -—' h
C = —,—6 JP .
2 2310

Similarly, (2.118) is inserted into (2.99) producing C3:

2

. . ism t+ sin t
03: EQerte—JMh-t) ‘1 q q q

221?

 

_ 2 '
7% sin qt + 2jpcos qt — qsin qt

sin qt + 2jpcos qt — qsin qt

(2.119)

(2.120)

(2.121)

The constants C4 and C5 can be found using a similar procedure. However, they are

not needed to solve the canonical problem and thus are left undetermined. Therefore

in summary,

0, Z flaeiph+fme—2ip(h—t)eiph(m

 

2119 2119

[#0 —' h
C = —,—8 Jp
2 29p

06 ' t — '- h—t
C = —.——er6 M )R
3 gm ( )
C6 = [#0

_ 2
7% sinqt + 2jpcosqt — qsinqt

where

_ 2
—p—sin qt + qsin qt
R: ‘1

 

_ 2 '
—qE- sinqt + 2jpcos qt — qsin qt

24

(2.122)
(2.123)

(2.124)
(2.125)

(2.126)

(2.127)

2.1.6 Electric Field Derivation

With the constants now known, the fields on each side of the dielectric slab may be
found. For the region above the slab, z 2 t, the solutions for the potential in regions

1 and 2 are combined as follows:

Ag;(ky,z) : fie JPIZ hI + $er JP(2+h)e2JPtR. (2.128)
21p 2119

Substituting (2.132) into the inverse Fourier transform wave equation, (2.42), gives

00
1 1 _.-. _
—OO

00
i fle—jp(2+h)€2jptRejky1/dky, (2.129)
27r 23p

~00

OI'

 

—OO

00 . .

I , —]p|z—h| . —j[)(2+h—2t) ’

_’10 [_e efl‘yy + e p Rejkyy dky. (2.130)
p '

The electric field above the slab, z 2 t, is then determined from (2.60) yielding

 

I Doejkyy
152(22): 4:0” / —p—[e‘JPIz‘hI+e"JP(Z+h‘2t)R]dky. (2.131)
—00

When 2 < 0 the vector potential is determined by (2.66). Substituting (2.125)

into (2.66) gives

- 1 21)»? ,—jI)(h-t)
4208. z) = 2 “06 e . (2132)
—Z;— sinqt + 2jpcos qt — qsinqt

 

25

OI‘ .
- [HOeJMz-hfi)

 

 

 

 

 

—%— sin qt + 2jpcos qt — qsin qt
Implementing the inverse Fouier transform, (2.42), results in
I I ’jp(z—h+t) jk'yy
42(9, 2:) = 7:9 2 ”0‘3 e dky. (2.134)
—00 —% sin qt + 2jpcos qt — qsin qt
Applying (2.60), the electric field for z < 0 is
00 . -
_ i, 1 1 Jp(z—h+t) Jkyy
E2(y. z) = —]—“2"7rfl 2 “06 e (12,). (2.135)
—00 ~%— sinqt + 2jpcos qt —- qsinqt
In summary, for the canonical problem seen in Figure 2.2, the electric: fields is found
in the regions above and below the dielectric slab:
00 .
—I , , ejk-yy _ .' _
Emu/,2) : “Oi-U / _[e JplZ h] +
477 p
—00
(WWW—20R] city, 2 2 t (2.136)
_ ' 1 ”(Z-[1+0 kyyj
Ex(y, z) : i2? 2 e 8 dry 2 < 0. (2.137)
—00 —%— sin qt + 2jpcos qt - qsin qt

2.2 Scattered Field from Conducting Strip

With the impressed fields on either side of the dielectric slab determined, the current.
across the infinitesimally thin conducting strip can be found. From the current, the

scattered field can be derived at a position above the slab.

26

2.2.1 Integral Equation Derivation

An integral equation is determined to obtain the current on the conducting strip [10].
Along the surface of the conducting strip the tangential electric field is equal to zero.

Given that the tangential component of the electric field is zero, from Figure 2.3 it is

Ej(y,z) E§(y,z)

 

 

> z=4
-w P w
Figure 2.3: Conducting strip parameters.
determined
1320.2) = —E2<y,z>. (2.138)

The strip is assumed to be infinitesimally thin, therefore the currents on the top
and bottom are nearly adjacent, which allows them to be summed together. The
integral representation of the scattered electric field in terms of the two-dimensional

Green’s function is

Eflyiz) = f / 01(9’,z’ly,z)J2(y’)dy’dz’. (2.139)

The summed current, J$(y), is independent of the :1: direction, which makes it possible

to write the surface current as a volume current density:

J$(y, Z) 2 J31; (5(2 - 25) (2.140)

27

Substituting (2.140) into (2.139) yields

'11)
/ Guy/.2311}, Zleszll/Wy, = way, a). (2.141)

y ——w

The integral equation in (2.141) can be solved using the method of moments
(MOM). This approach involves expanding the unknown current as a linear combina-

tion of basis functions:

N
J3$(y) = Z anfn(y). (2.142)
n=1

Substituting (2.142) into (2.141) gives
N “1’
2: an / 01(y’,zsly,zs>fn<y’>dy’ = —E§-(y,zs>. (2.143)
n=1 —w

For simplicity the basis functions fn are chosen to be pulse functions:

lryn-%Sy’S3/n+%

fn(y) = P7101) = (2.144)
0 : elsewhere.
In (2.144),
3177. = —w + (n - g) A, (2.145)
and
A _ 22 (2 146)
— N o u

The quantity A is the width of the partitions, N is the total number of partitions
used, and yn is the center of a partition.
The method of collocation (or point matching) is used to convert the integral

equation into a system of linear equations. By point matching at y = ym, a set of

28

linear equations are obtained:

A
Z (In / Gl(y,,Zsly-7n,23)dy, = —E:tr(y7n,23), m =1,2,...,N, (2.147)
n=1
3171—79—
01'
N
Z anAynn = bm. (2.148)
n=1

The variable an in the matrix equation (2.148), represents the current on the con-
dUcting strip.

When a field impinges on an object, current is induced along the surface. Scattered
fields are then produced from this induced current. For this canonical problem, the
incident field on the conducting strip is the transmitted field from an electric line
source above a lossy dielectric slab. Therefore giving,

00 . - ,
b —jw]/t0 ejp(z—h+t)le~y?/m
m = _—

 

dky. (2.149)

2 2 . . .
—oo —% sm qt + 2jpcos qt — qs1n qt

The scattered field from the strip is actually the impressed field produced by an

electric line source in the presence of a lossy slab. Thus yielding

 

 

oo . I
— thyUJm-y) . /
. U.) 6 _ _
61(y’.z'|ym,2)= :0 f [e 3142 2|+

7‘ 1)

—OO

. ~ ,—
ReJP(~+Z 2t)]dky. (2.150)

The source point in the equations for the impressed and transmitted fields derived in
section 2.1.6 is in reference to the front of the dielectric slab. Therefore the source
point in (2.149) and (2.150) also have to be referenced from the front of the slab.

This results in z, = t — 2:3. Since (2.148) is defined for all points on the surface of the

29

conductor, the observation point also results in z = t — zs. Rewriting (2.147) yields

 

 

 

A
971+? 00 . . /
‘ltow BJkyyme—Jkyy 7.1)27
Amn = —— [1+ Re- 4.9] dkydy. (2.151)
47r p
A —00
Lyn—'2'
Rearranging (2.151) gives
00 - '
_ _ Jky'ym .
AW, 2 "0‘” 6 [1+ ReJP22‘8] Xdky, (2.152)
47r p
—00
where
yn+%' I
X = / e‘Jkyy dy’. (2.153)
yn—%
Performing the integration in (2.153) yields
A
e_jkyyl girl-‘2-
X = ——7— . (2.154)
3 ‘31}
3171—7123—

Inserting the limits of integration for the variable y, gives,

. A . - . A
_e—Jkyynejky?’ + e‘JkyyneJk'y'Q‘
X = , . (2.155)
Jky

 

Simplifying further gives

'. A _,-, A
—jr~. me”? -e ””7
X = e y , , (2.156)
3163/

30

which can be reduced to

. . - A
X = e—jkyynm,
Jky
or A
X = e_jkyynASlnky—2—
k A '
1’17

Substituting (2.158) into (2.152) produces

 

 

oo . .
_ _Jk'l 3m J}? 31771 . -
Amn = [10w / e J 8 y [1 + ReijzS]

p A
-00 km

2.2.2 Amn Evaluation

The Amn integral actually consist of two separate integrals:
Amn = Umn + an,

where 00 /’
wow / e-Jkyyneykyym Asin kyTAdey
A i 9
471’ p ky‘g’

—OO

Umn '—

 

 

and

 

 

oo - -
_ A: ., . .. . . A
_”0(, f e J yyneflxyymAsmky? ejPQZs

477 p

(1101 .
kg?— J
—00

an =

A sin ky%- dk

y.

(2.157)

(2.158)

(2.159)

(2.160)

(2.161)

(2.162)

It is computationally expensive to integrate a function from negative infinity to posi-

tive infinity. Since portions of the integrand are even about it’s axis, the integrations

of Umn and an are simplified. First employing Euler’s relation to (2.161) results

31

in

 

 

oo _ A
—/1,0wA cos ky(m - n)A S111 [Cg/72-
Umn = +

4 A
7r —00 p [Cy-2—
" - k A

J sin 1.1/(m — n)A Sln 31-?-

p 16,179

 

dry. (2.163)

Given the definition of a sine function, (2.163) is rewritten as

 

 

UNIT), =

00
—/10wA / cos 163/(m — n)A

'ck A+
sm " —
471' p y2

0
j sin ky(m — n)A
p

 

A
sinc ky—2—dky. (2.164)

In (2.164), the sine and 19 functions are even about the vertical axis. When these
functions are multiplited by a cos this results in an even integrand. While mulitplying

by a sin function results in an odd integrand. Therefore this yield

 

 

 

 

 

 

oo
— k( — )A A
Um”: #201511 / COS 9(1) m n) alleging, (2.165)
0
or
A 00 k( m n) )A sin kyA
— 1 52 cos —
Umn— _ ’ 0 / y( A7111.) (2.166)
0 k9?
Through the same reasoning (2.162) results in
00 A
— ,l;( —- )A sin k
an= “20cm / COS '(y m n) :7 ReJP2~Sdly . (2.167)
0 kg?

32

2.2.2.1 Umn Evaluation
The integrand in (2.166) has the form of

OO

sin (11: cos 6.7:
Umn -—— ——-—(f$.

0 :rVk2—562

(2.168)

When a: = k the denominator of the integrand goes to zero, which is a singularity of

the function. Therefore a numerical method [12] is considered which splits Umn up

into two parts, thus removing the singularity:

k f( > 00 f( >

11} CE

U7nn=/————(1I+/———d$,
0 k2_$2 k ‘/k2_$2

where
sin arr cos b1:

f0?) —

13

Rewriting (2.169) gives,
U77"), = [it + 13’.

Looking at I i‘, the square root in the denominator can be factored giving

Vk+$

k _fleL
Igzo/(T/Efildx,

OI'

 

Now, use the substitution

(2.169)

(2.170)

(2.171)

(2.172)

(2.173)

(2.174)

01‘

t: x/k — x, (2.175)

with the differential

(1:1: = —2tdt. (2.176)

Inserting (2.175) and (2.176) into (2.173) and changing the limits of integration results

 

in
u 9(k — t2)
11 = —T—-(—2tdt) (2.177)
«E
or
fl
1? = 2 / g(k -— t2)dt. (2.178)
0
In the equation for I u, :1: is always larger than k. given the limits of integration.
Therefore,
00 f( >
I“ = / ‘E (1.1:, (2.179)
k —j‘ /:1:2 _ k2
or

15‘ = j (2.180)

00
/ —\/_f2—(_L——2d:r.
k a: — k
The sign on the square root function in the denominator of (2.179) is chosen because
this function represents the variable p in the original equation. In (2.68) the sign on
p is determined through physical reasoning.

The integrand in (2.180) oscillates substantially and only decays as 1/232; thus it

is difficult to integrate. Using Kummer’s method [13],

CX) OO
. f ,1 .
12‘ = 7 1/ [7,9375 — 16(4)] 4:: +1 [1144114, (2.181)

34

where fa(:r) is asymptotic expansion of —ié\/—'—T—)—k—2, the integration is easier to com-
m _—

pute. This is useful only if the oscillations in the integral cancel out at high multiples

of k, and the second term can be integrated in closed form. Further manipulations

of (2.181) yields,

(2.182)

 

0° ,G—_—2 2
u2j/f(l -\/l37:-_kf2fa(l )dr+j/fa(a:)dx
1.

Similarly to the simplification done for 11‘ , the square root function is split up:

ooff37)‘ V $2—k2fa(17) oo

 

 

 

15‘ = j/ H d1: +j [fa(:r)d:1:, (2.183)
16 .
or
00 ~ oo
23' / $111: + j / fa(:r)dx. (2.184)
16 k
Next, use the substitution
x=k+k? (2mm
t = a: — k, (2.186)
with the differential
d3: = 2tdt. (2.187)

By substituting (2.186) and (2.187) into (2.184) gives

- _, t
15‘ = j / LP2tdt+ j / fa(;1:)d;1~., (2.188)
16 is

01'

OO 00
15‘ =2j /g(k+12)dt+j/fa (:1:)d;1:.
k k

In a general form,

00
1% = 2j/§(k+t2)dt+jIA,
k

where

00
IA = [fa(:1:)d:c
16

(2.189)

(2.190)

(2.191)

With the equation for 1%“ in manageable form the the asymptotic expansion fa(;1:)

is chosen to achieve the requirements for the Kummer’s method. By choosing,

sin (1:1: cos bx

fa(~’13)= fi2+$2

,

and 13 to be arbitrary, I A can be rewritten as

M=%-k,
where
00 . b
s1na:1: cos :1:
I = / —————d:r,
B 32 + $2
0
and

I /sina:1:____2_cosb:1:d$
0: fl2 + 9:2

The equation for [C can be computed numerically. However I B is

(2.192)

(2.193)

(2.194)

(2.195)

more difficult.

To calculateIB, exponential integrals are used. Writing I B in terms of exponential

integrals yields

[B = 5 [E [(a — 0131+ E [(a + 1113]],

36

(2.196)

where

E(;zt) = e_xE.,;(1:) — emEi(—:1;) (2.197)
and :1:
—-t
EZ-(x) = / fit—111. (2.198)

2.2.2.2 an Evaluation

When considering an in (2.167), the upper integral limit is changed to a constant,
A. This constant is chosen large enough in value so that the integrand has decayed

close to zero and no longer significantly contributes to the integral. This results in

A
k
0

.2 _ 2
’"0 kg
The constant A must be larger than k0, since it is at kg = 160 that the integrand
blows up and after where the integrand decays to close to zero. Therefore A, is chosen
to be a multiplicative factor of 160.
A numerical method [12] is used to work around the singularity when kg = 1:0.

Splitting the integral at k0, (2.199) becomes

ll0
k k
k2 _ k2 k2 _ k2
0 0 y [co 0 2U
01‘

37

Evaluating I i) by first factoring the square root in the denominator gives,

k0 [M]

1v=/—__
1 0 meg—kg

dry,

01‘

’10
Iv [9 (1%)]

”Om

for simplicity. Then by using substitution

city,

A _. 2
ky—ko—t

OI‘
t: (LIED—kg,

along with the differential
dky = —2tdt,

in (2.203) produces

Simplifying further results in

m
11) = 2 / g(k0 — t2)dt.
0

(2.202)

(2.203)

(2.204)

(2.205)

(2.206)

(2.207)

(2.208)

Knowing that the constant A is larger than 360, the second part of an can be

38

rewritten as

 

, k
I?“ _ f( 2”) 21111,).
‘3 ky — k0

Then, by factorng the square root in the denominator produces,

f‘

A [Jill/L]

1; —1.~

Using substitution
or

along with the differential
dky = 2tdt,
11 (2.210) result in

()(k t2)
—0—(t+ ——(2tdt).

Then with further simplification,

A—ko

[g = 2j / 9050 + t2)dt.
0

39

(2.209)

(2.210)

(2.211)

(2.212)

(2.213)

(2.214)

(2.215)

2.2.3 Matrix Equation

With Amn put into manageable form and

—jw1#0 00 ejp(z—h+t)ejkyym
bm = T

 

2 dky, (2.216)
——oo —2q— sin qt + 2jpcosqt — qsin qt

the current in the matrix equation of (2.148) can be derived. Thus the current at
the center of each of the partition along the strip, an, is computed. These current
values then can be treated as individual current sources with a width equal to that of
the partition width of the conducting strip. The fields produced from an individual
current source in the presence of a dielectric slab was derived in (2.137). Therefore
the scattered field from strip can be computed by summing up the fields produced by
each of the individual line sources. With this derivation the reflected and transmitted

electric fields for the canonical problem seen in Figure 2.1 is solved for.

40

Chapter 3

Numerical Results

Computation of the fields scattered by the strip is carried out using the FORTRAN
programming language. The position of the line source is set and represents the
location of a transmitter. The scattered field is computed at an adjacent point,
representing the position of a receiver. As seen in Figure 3.1, the transmit / receive
pairs. are moved parallel to the barrier and the fields computed from 20 MHz to 16

GHz at 20 MHz increments for each pair.

 

21 different transmitter/receiver positions _
I Equally spaced over 5 m |Z-0.5m

...XX XX XX XX XX XX XX XX XX XX XX XX XX XX XX...

 

 

 

 

 

 

80 z z=o.1524m
)5, =6 0=.001$/m g
V =
80 y 2 0
, .15 m 1 Z=-O.8m

Figure 3.1: Position of source/ observation pairs.

The total field is computed at 21 different transmitter/ receiver locations in the

frequency domain. At each of these positions, the frequency domain responses are

41

transformed into the time domain using the inverse Fourier transform. The time
domain data fields from all 21 transmit/ receive positions are used to construct an
image of the target. An observation region is created using a 2-dimensional grid of
points, and the propagation time for a wave propagating from each of the 21 different
transmitter/receiver positions to that point is calculated. Given these propagation
times, the intensity values from the time domain data are summed together creating

an image of the target.

3.1 Computation of Incident and Scattered Fields

The permittivity, conductance, and thickness of the slab, as shown in Figure 3.1, are
chosen to those of a concrete wall [14]. Multiple reflections from the conducting strip
to the back of the dielectric slab will produce ghost images of the target. Therefore to
decrease the merging of the ghosts and target responses, the strip is located sufficiently
distant (z = —0.8 m) from the back of the dielectric slab. The width of the conducting
strip is chosen to be relatively small at a width of 0.15 m., to test the ability of the
image technique to resolve small targets.

The conducting strip is partitioned into narrower strips, with each of these strips
having a constant current. The MoM technique is used to find the scattered field from
the conducting strip. The number of partitions on the strip is set to 20 per wavelength,
and is thus dependent on the frequency used. This number can be computed using
the formula

2111f

N = 20—— 3.1
C < >

where 211) is the width of the strip, f is the frequency and c is the speed of light. The
minimum number of partitions allowed is 10, regardless of frequency.

The integrals needed to find the incident and scattered fields are computed using

the FORTRAN subroutine DQAG, which is part of the QUADPACK computational

42

package [15]. The computational package includes approximations to integrals of the

form

b
/f(:1:)w(:1:)d:1:. (3.2)
a

The weight function w is used to incorporate known algebraic or logarithmic singu-
larities, oscillations, or to indicate that a Cauchy principal value is desired. DQAG
is a general purpose integration routine with a high efficiency and w(:1:) = 1 in (3.2).

The routine subdivides the interval [a, b] and uses a (216+ 1) -point Gauss-Kronrod
rule to estimate the integral over that range [11]. The Gauss-Kronrod rule is an
adaptive Gaussian quadrature method in which, for numerical integral routines, the
error is estimated based on the evaluation at special points called Kronrod points.
With this suitable point selection, the x-coordinate from previous iterations can be
reused as data for the new set of points. A Gaussian quadrature would require a new
x-coordinate at each of the iterations; this is important when a specified degree of
accuracy is required but the number of points needed to achieve this accuracy is not
known in advance.

The calling sequence for DQAG includes the function to integrate F, the lower
limit a and upper limit b. The user defined absolute error e is ERRABS, and the
relative error ,0 is ERRREL. The DQAG routine returns two numbers of interest:
RESULT and ERREST, which are the approximate integral R and the error estimate

E, respectively. The relation between these numbers are as follows:

b b
/f(:1:)d:1:—R _<_ESmax e,p /f(;1:)d:1: (3.3)

a a

0‘200 and the relative accuracy is

In these simulations, the absolute accuracy is 1
chosen as 18 decimal places. Since the limits of integration are infinite, it is necessary

to find a point to truncate the integrals that provides a desired accuracy, hence

43

determining a and b. Since the integrands of (2.136) and (2.137) are even about
ky, the lower limit of the integral can be set to 0 and the result multiplied by 2.
The integrand decays rapidly once ky > k0. When ky > 160, values of p become
imaginary in (2.50) producing an exponential decay in the integrand. Therefore, a
good means. for determining the integration truncation point is a multiple of 160. In
these simulations the multiple is chosen to be 10.

The scattered field from the conducting strip is found by first solving the matrix
equation (2.148). As described previously, Amn is separated into two integrals. Using
the DQAG integral routine, bmn and an are computed for the specified degree of
accuracy. In addition, Umn requires two more subroutines to evaluate E1(x) and
E1(x), the exponential integrals found in the evaluation of Umn. The matrix equation
is of the simple form Ax = b, which is solved using another subroutine called SOLVEC
using the decomposed matrix from the routine DCOMP. The solution is the complex
current at each partition on the conducting strip. Finally, the field is determined at
the receiver position by computing the scattered field from the known current on the

conducting strip.

3.2 Imaging Analysis

3.2.1 Frequency to Time Domain Transformation

The field is computed at 21 positions above the slab between -2.5 m to 2.5 m with a
spacing of 0.25 m, as seen in Figure 3.1. At each of these locations, the incident and
scattered field is computed in the frequency domain between 20 MHz and 16 GHz,
at a step size of 20 MHz. The incident and scattered fields are computed separately,
and the frequency domain responses are transformed into the time domain using the
visual basic program WaveCalc (supplied by John Ross and Associates). An example

of the frequency domain field is shown in Figure 3.2 for the transmit/receive pair

44

centered above the conducting strip.

 

1000 -

800 '

 

 

 

 

 

 

 

 

 

Magnitude

 

 

 

 

 

 

 

0'1...1...1...lm..1...1...1...1...1
0 2 4 6 8 10 12 14 16

Frequency (GHz)

 

Figure 3.2: Magnitude of the frequency domain scattered field at a position centered

above the strip.

Computing the fields over a certain frequency range is equivalent to evaluating
the total fields over frequencies ranging from zero to infinity and truncating them
using a rectangular function. Therefore, when the data is transformed into the time
domain, the rectangular function is transformed into a sine function, producing an
abundance of oscillations on the edges. To eliminate some of these oscillations, the
data is windowed in the frequency domain using a Gaussian modulated cosine (GMC)
instead of a rectangular function. This allows for a steady roll off of the data, so that
there are no abrupt changes, which results in few oscillations in the time domain. For
these simulations, the GMC is centered at zero frequency with an equivalent pulse

width of 0.12 ns. Figure 3.3 shows the weighted frequency domain data.

45

 

 

Magnitude
8
. , . .

40;

20:

 

 

I..1.l.1...11..1...1...1...1...1
O 2 4 6 8 10 12 14 16

Frequency (GHz)

 

Figure 3.3: Magnitude of the frequency domain scattered field weighted with a GMC.

Finally, the weighted frequency domain data is zero padded and WaveCalc is used
to transform the data into the time domain. Zero padding adds additional points at
high frequencies and results in a finer sampling in the time domain. The final time

domain response is shown in Figure 3.4.

46

 

81:1011 :'

6x10" :-

'TT

4x1011 _

N
X
.1.
‘3.
I

 

Magnitude
O
l
_aL.
.1:—
1r
4-

 

I l I

-2x10"

-4x10‘1

TTITI

 

III

51:1011

 

 

'1....1....1..411....1....1
0 10 20 30 40 50

Time (ns)

 

Figure 3.4: Time-domain scattered field at a point centered above the conductive

strip.

The first event seen in Figure 3.4 represents the first wave interaction with the
conducting strip. The following spikes are multiple reflections between the dielectric
slab and the strip. Note that Figure 3.4 does not include the incident field. This
field, shown in Figure 3.5, is computed separately and shows much stronger spikes

representing the field reflected by the slab (the flash).

47

 

2x1013 :'

1111013}

Magnitude
53.
i"

-2x1013 _"

 

T'I‘I

 

-31<1013

 

 

l l l j l l l l I A l L A A; L 14 1 1 1 l l 1 1 #1

o 10 20 30 4o 50
TIme (ns)

 

Figure 3.5: Time-domain incident field at a point centered above the conducting strip.

3.2.2 2D Image Construction

The time domain data fields from all 21 transmit / receive positions are used to con-
struct an image of the target. An observation region is created using a 2-dimensional
grid of points, and the propagation time for a wave propagating from each of the
21 different transmitter/receiver positions to that point is calculated. Given this
distance between the transmit/ receive position and the grid point, the approximate
propagation time is

T = —\/€7. (3.4)

C

Where 7' is the travel time, d is the distance from the transmitter / receiver position to
a particular point on the grid, 6 is the speed of light and 67- is the relative permittivity
of the material through which the wave travels. The actual travel time is somewhat

difficult to compute since a ray entering the lossy slab undergoes refraction, and thus

48

doesnt follow the straight line path between the transmit / receive point and the grid
point. However, since the slab is thin, the straight line path gives very accurate
results. For the portions of the path in free space, er is unity, while for the portions
inside the slab ET is taken to be the real part of the complex permittivity of the slab:
CT = 6.

Using the computed travel time, the value of the total (incident plus scattered)
field is found from each time domain waveform, and the results for all 21 transmit / re-
ceive positions are added to determine the image intensity at this particular grid
point. The image intensity is highest at points corresponding to places of strong
reflection, such as the front of the slab and the target. Thus, strong image intensity
indicates a scatterer. Sometimes this is called a scattering center approach since it
localizes the points on the scatterer that produce the largest reflections.

Figures 3.6 and 3.7 display the two-dimensional intensity image within a box
that includes a portion of the dielectric slab and the conducting strip, using the
time-domain data from 21 transmit/ receive pairs. The axes of these two figures are

distance in meters.

49

H~AL11 1'11 .14‘41‘4- cc .1 "J1" ll.
“CVW 1 - . .. ..
(_p § \ - s

’ '- 4' 5A5

 

Figure 3.7: Proportionally correct mage of lossy dielectric slab and conducting strip

with scale x1013.

50

Figure 3.6 shows a compressed aspect ratio resulting in a square figure, while
Figure 3.6 shows the proper aspect ratio. As seen in Figure 3.1, the back of the lossy
dielectric slab is positioned at z=0 m, while the front is at 220.1524 m and the strip is
at z=-0.8 m. Comparing this figure to Figure 3.6 visually shows the correct positions
of the objects analyzed.

Figure 3.6 shows how the image of the target and the slab are built up from the
time responses at each of the transmit / receive pairs. A strong peak in a time response
produces an arc in the image. These arcs overlap at points where a scattering center
exists, since the reflection from the scattering center appears in all of the observations.
Thus, the brightest spots in the image correspond to the strongest scatterers. In

particular, the conducting strip target can be clearly seen.

3.2.3 Image Quantification

The sharpness of the target image depends on the parameters of the simulated radar
system, including sampling rate, pulse width, etc. To determine the effects on the
image of altering certain parameters, a way of quantifying the sharpness of the con-
ducting strip image is needed. To do this, an image radius is defined, which is the
radius of an equivalent circle that has an area equal to the region that contains the
majority of the image intensity. An observation box is centered on the conducting

strip as shown in Figure 3.8.

51

 

-2 -1 0 1 2

Figure 3.8: Area in which sharpness of conducting strip is determined with scale

x1013.

The size of the box is arbitrary, but it should completely contain the strip and
should not be too large so that other artifacts in the image are excluded. In this
research, a box size of 0.3 by 0.46 111 is chosen. Inside this box the position of the
center of pixel intensity is determined. This approach is much like the formula for

finding the center of mass of an object. The coordinates of the center are given by

”y
g: fiyi
yo = $ (3.5)
,2 f4
1:1
712
Z fizi
20 = _1: , (3.6)

In these equations y and z are the pixel positions, 71y and 712 are the number of pixels
in the y and z direction respectively, and f represents that particular pixel intensity.
The radius of intensity, R, is found by considering only pixels with intensities

above 0.5 times the maximum intensities of all pixels within the box, which are used

52

in the formula 71

12 f1 (91-1/0)2+(21-20)2
R: _

 

n (3.7)
2 fz'
z=1

where n = nynz. The image radius is the radius of a circle that contains the higher

pixel values and gives a single quantity through which the sharpness of the strip image
may be described. A large image radius corresponds to a blurry image and a small
radius to a sharper image. An example of an image radius is shown in Figure 3.9 for

a relatively large pulse width of 2.0 ns.

0.5

 

-2 -1 0 1 2

Figure 3.9: Circle of radius equal to the image radius centered around the conducting

strip with scale x10“.

53

Chapter 4

Parameter Analysis

A parameter analysis is used to compare the performance of the ultra—wide band
short-pulse and stepped frequency radar systems. Simulated data is computed at
21 different transmitter/ receiver positions, images are created for different parameter
sets and the sharpness of the images is analyzed. When the pulse width, sampling
rate, SNR, or digitization of the sample data is altered in both the frequency and
time domain, the image radius changes allowing for optimization of these parameters

in an actual prototype.

4.1 Time Domain Analysis

The significant parameters that are investigated during simulations of the time-
domain system are the incident pulse width, sample rate, SNR, and digitization.

Each of these is considered here in turn.

4.1.1 Pulse Width

The pulse width in the time domain is the duration of the pulse emitted from the

transmitter of a UWB radar system. Longer durations of the pulse width cause the

54

image of the conducting strip to become blurred. Additionally, during transmitting
the receiver is non-active to prevent cross coupling of the transmitted pulse. If the
duration of the pulse is longer than it takes for the start of the pulse to reflect from an
object and return to the receiver, this signal is not received, producing false results. To
simulate this type of parameter change, the width of the Gaussian modulated cosine
window function is altered in the frequency domain. With the geometry parameters

of Figure 3.1, values of the computed image radius are shown in Figure (4.1).

 

1-l —L .5
O N «h
IIUII'IIIIII

Image Radius (cm)
0) on
' I ' ' I

 

 

 

4:-

l
2—

'1 1 1 1 1 l 1 1 1 14 1 1 1 1 I 1 1 1 1 l 1 1 1 1 [m1 1 1 l
0 0.5 1 1.5 2 2.5 3

Pulse Width (ns)

Figure 4.1: Effect of pulse width on image radius.

This figure demonstrates that increasing the pulse width enlarges the image radius
which corresponds to an increase in the blurring of the strip. Between 0.06 ns and 0.10
us the image radius unexpectedly decreases; this probably is because the threshold
value of the pixel intensity used to compute the radius is set too low.

It is helpful to visually associate the reduction of image quality as the pulse width

increases with the increase of image radius. The images are displayed in Figure 4.2—

55

4.11. When the pulse width reaches 1.5 ns, Figure 4.8, the conducting strip nearly
fades into the background. In Figure 4.11, where the pulse width is 3.0 ns, the image
of the conducting strip appears to have moved away from the actual position of the
strip. The size, position, intensity, and shape of the strip can not be determined when

the pulse width is too large.

0.5

-1

 

-2 -1 0 1 2

Figure 4.2: Conducting strip image constructed using a pulse width of 7‘ = 0.06 ns.

56

 

-2 -1 0 l 2

Figure 4.3: Conducting strip image constructed using a pulse width of 7' = 0.10 ns.

57

 

-2 -l 0 l 2

Figure 4.4: Conducting strip image constructed using a pulse width of 7 = 0.16 ns.

58

 

-2 -l 0 l 2

Figure 4.5: Conducting strip image constructed using a pulse width of T = 0.26 ns.

59

 

-2 -1 0 l 2

Figure 4.6: Conducting strip image constructed using a pulse width of T = 0.50 ns.

60

 

-2 -l 0 1 2

Figure 4.7: Conducting strip image constructed using a pulse width of 7' = 1.00 ns.

61

 

-2 -1 0 1 2

Figure 4.8: Conducting strip image constructed using a. pulse width of T = 1.50 ns.

62

 

-2 -1 0 l 2

Figure 4.9: Conducting strip image constructed using a pulse width of 7' = 2.00 ns.

63

 

-2 -l 0 1 2

Figure 4.10: Conducting strip image constructed using a pulse width of T = 2.50 ns.

64

 

-2 -l 0 1 2

Figure 4.11: Conducting strip image constructed using a pulse width of T = 3.00 ns.

4.1.2 Sampling Rate

The sample rate in the time domain represents the number of samples taken per
second from a continuous signal to form a discrete signal (also described in terms
of the sampling period, which is the time between samples). The simulated data is
computed at a high sampling rate, and the affect of lower sampling rates is easily
determined by discarding data points. Figure 4.12 shows the dependence of image
radius on sampling rate. As the number of samples per second decreases, the image

radius increases and the conducting strip becomes blurrier.

65

 

‘
1h
It'llrll

_L
N

.1;
O
' 1

Image Radius (cm)
a: 00
. , . . ,
El

 

 

l 14 l l I I I l l L lJ 1 I1 I I l1 1 1 L1 1 l I l ngl

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
Sampling Period (ns)

 

Figure 4.12: Effect of sampling period on image radius.

As seen with the affect of pulse width, the smallest value of sampling period does
not produce the smallest image radius. This is again a result of the difficulty in
choosing an appropriate threshold for selecting which pixel intensities are to be used
in computing the radius. For visual inspection, Figures 4.13 - 4.22 show the images
for each of the sample rates shown in Figure 4.12. As seen in Figure 4.17, when the
sample period reaches 0.40 ns the image of the conducting strip is so unclear that an
accurate detection cannot be produced.

In addition to blurriness, another concern is the existence of ghost images of the
conducting strip. With the time-domain system, ghosting is due to multiple reflections
within the barrier or between the target and the barrier, which produce additional
peaks in the scattered field isignal. These peaks add when constructing the image,

producing ghost images at regular intervals behind the actual image. In Figure 4.17

66

it is seen that for a large sample period (0.40 ns, which is a sampling rate of 2.5
GHz), the intensity of the ghost image is larger than the actual conducting strip,
which would result in an incorrect location of the object behind the slab. When the
sampling period reaches 1.60 ns, Figure 4.22, the conducting strip is not visible and

no analysis is possible.

0.5

-0.5

-1

 

-2.5

-2 -1 0 1 2

Figure 4.13: Conducting strip image constructed using a sampling period of At =

0.05 ns.

67

0.5

 

-1.5

-2 -1 0 1 2

Figure 4.14: Conducting strip image constructed using a sampling period of At =

0.10 ns.

68

 

-2 -l 0 1 2

Figure 4.15: Conducting strip image constructed using a sampling period of At =

0.20 ns.

69

0.5

-0.5

-1

 

-1.5

-2 -1 0 1 2

Figure 4.16: Conducting strip image constructed using a sampling period of At =

0.30 ns.

70

0.5

 

-l.5

-2 -1 0 1 2

Figure 4.17: Conducting strip image constructed using a sampling period of At =

0.40 ns.

71

 

Figure 4.18: Conducting strip image constructed using a sampling period of At =

0.50 ns.

72

 

Figure 4.19: Conducting strip image constructed using a sampling period of At =

0.60 ns.

73

0.5

-1

 

-1.5

-2 -1 0 1 2

Figure 4.20: Conducting strip image constructed using a sampling period of At =

0.70 ns.

74

0.5

 

-1.5

-2 -1 0 1 2

Figure 4.21: Conducting strip image constructed using a sampling period of At =

0.80 ns.

75

0.5

-1

 

-1.5

-2. -1 o 1 2

Figure 4.22: Conducting strip image constructed using a sampling period of At =
1.60 ns.

4.1.3 Signal-to—Noise Ratio

An essential consideration in any radar system is the signal—to—noise ratio (SNR)
appearing at the receiver. As the SNR is reduced, it is anticipated that the target
image will become less distinct. To examine this effect, random Gaussian noise is
added to the simulated time-domain data using the program WavCalc. The noise
level is relative to the largest peak in the waveform (which is the flash from the
barrier), and thus the SNR values described here represent the decibels (dB) down
from the peak value of the time signal. The relationship of image radius to SN R is

shown in Figure 4.23.

76

 

 

 

 

D

8~El

E.

3

(06'

3..

:6 C]

(U

m.

m—

CD4

CU”

E

— C] D

' Cl

2—

. U D D
. Cl
1‘..|..ilii.li..l.i.l.rnl...I..i|.n.l
2 4 6 8 1O 12 14 16 18 20

SNR(dB)

Figure 4.23: Effect of noise on image radius.

While the sharpness of the image is dependent on SN R, it is seen that the noise
must be around SNR = 10 dB before image degradation occurs. This is because
reflection from the target is quite large, and the time-domain peak extends above the
noise until the noise level becomes as large as the reflection.

The target images as a function of SNR are shown in Figures 4.24-4.33. Though
the image radius at SNR = 4 dB is not as large compared to the radii produced by
high sample periods and pulse widths, at the low SNR the conducting strip still blends
into the background, as seen in Figure 4.32. As the noise increases, the background

becomes uniform and no images can be discerned.

77

 

-2 -1 0 1 2

Figure 4.24: Conducting strip image constructed using a SN R = 20 dB.

78

 

-2 -1 0 1 2

Figure 4.25: Conducting strip image constructed using 3. SNR = 18 dB.

79

 

-2 -1 0 1 2

Figure 4.26: Conducting strip image constructed using a SNR = 16 dB.

80

 

-2 -1 0 1 2

Figure 4.27: Conducting strip image constructed using a SNR = 14 dB.

81

 

-2 -1 0 1 2

Figure 4.28: Conducting strip image constructed using 3 SNR = 12 dB.

82

 

-2 -1 0 1 2

Figure 4.29: Conducting strip image constructed using a SNR = 10 dB.

83

 

-2 -1 U 1 2

Figure 4.30: Conducting strip image constructed using a SNR = 8 dB.

84

 

-2 -1 0 1 2

Figure 4.31: Conducting strip image constructed using a SNR = 6 dB.

85

 

-2 -1 O 1 2

Figure 4.32: Conducting strip image constructed using a SNR = 4 dB.

86

 

-2 -1 I] 1 2

Figure 4.33: Conducting strip image constructed using 3. SN R = 2 dB.

4.1.4 Digitization

In a time-domain radar system a fast A/D converter is used to sample the time
signal, digitizing the sample using a certain number of bits. Using a low number of
bits results in losing smaller signals by rounding, and poor accuracy in larger signals.
It is anticipated that dependence of the sharpness of the image on the number of bits
used in the A/ D converter will be fairly consistent until the number of bits becomes
too small. At that point the image should degrade rapidly.

In order to implement the digitization process, the simulated time data is dis-

87

cretized into bins of size (1 based on the number of A / D bits, according to the formula

R

d=2N_1. (an

 

Here R represents the range of the intensity and N the number of bits. Every point
in the time data must be put into a bin. This is done using the equation

bin 2 31min + round [W] d (4.2)

With this digitized version of the time data, images are created using different num-
bers of bits. The images are shown in Figures 4.35-4.44. Figure 4.34 shows how the

image radius depends on the number of bits used.

 

d —L d
O N -h
II II II

C)
l

I I

Image Radius (cm)
on
I

A
I I l I

 

 

 

12 1o 8 6 4
Bits

Figure 4.34: Effect of digitization on image radius.

As seen in the figure there are two jumps in image radius. From 12 to 9 bits

88

the constant image radius results from the conducting strip peak intensity values not
being effected by the rounding into bins. The jumps between 9 and 8 bits, and 7
and 6 bits result from a threshold bit value where the peak values new round poorly
into large bin sizes. From 7 to 2 bits the image radius is so high that decreasing the
bits can not decrease the image quality anymore. The subtle changes in intensity
become lost and soon, as seen in Figure 4.44, the strip completely disappears into
the background. This results from having too few bits to differentiate between the
strip intensity and the background. However, it is interesting to note that the ghost

images of the strip are actually removed from the image with decreasing bit numbers

since these weaker artifacts fall into the first bin.

0.5

 

-2 -1 0 1

2

Figure 4.35: Conducting strip image constructed using the number of bits = 12

89

 

-2 -1 0 1 2

Figure 4.36: Conducting strip image constructed using the number of bits = 11.

90

 

-2 -1 0 1 2

Figure 4.37: Conducting strip image constructed using the number of bits : 10.

91

 

Figure 4.38: Conducting strip image constructed using the number of bits = 9.

92

 

-2 -1 D 1 2

Figure 4.39: Conducting strip image constructed using the number of bits = 8.

93

 

-2 -1 U 1 2

Figure 4.40: Conducting strip image constructed using the number of bits = 7.

94

 

-2 -1 0 1 2

Figure 4.41: Conducting strip image constructed using the number of bits = 6.

95

 

-2 -1 U 1 2

Figure 4.42: Conducting strip image constructed using the number of bits = 5.

96

 

-2 -1 U 1 2

Figure 4.43: Conducting strip image constructed using the number of bits = 4.

97

 

-2 -1 El 1 2

Figure 4.44: Conducting strip image constructed using the number of bits 2 3.

4.2 Frequency Domain Analysis

The important parameters concerning the frequency domain radar system are band—
width (equivalent to pulse width in the time domain system), sample rate, SNR, and
digitization. The effect of these parameters on image quality is investigated by ma-
nipulating the simulated frequency domain data before the inverse Fourier transform

is taken.

98

4.2. 1 Bandwidth

Increasing the bandwidth of a stepped-frequency radar system increases the reso-
lution of the images generated by the system. To investigate this effect with the
simulated data, the frequency domain data is multiplied by a Gaussian modulated
cosine windowing function centered on 8 GHz. The equivalent frequency pulse width,
T , is changed to determine the bandwidth, (2, as given by

Q = —1—4ln 2. (4.3)

7”"

The dependence of image radius on bandwidth is shown in Figure 4.45.

 

 

 

 

10 l- U
f
_ El
- C]
A 8 -
E - a
U) E]
.2 D
a e -
a: C]
a . a
g s
_ 4 '- D
2 'Ci I . . . l . i . J . . . l . . . 4 a
100 80 60 4O 20

Bandwidth (%)

Figure 4.45: Effect of digitization on image radius.

The following figures show the decrease in image quality as the bandwidth in-
creases. Following the progression though each of the images, the conducting strip

and the lossy dielectric slab become progressively unclear. When the bandwidth is

99

at 5.81%, as seen in Figure 4.55, the conducting strip is nearly indiscernible from the
background; this makes locating the object behind the lossy dielectric slab almost

impossible.

0.5

-1

 

-2 -1 0 1 2

Figure 4.46: Conducting strip image constructed using a bandwidth = 110.38%.

100

 

-2 -1 U 1 2

Figure 4.47: Conducting strip image constructed using a bandwidth 2 36.75%.

101

 

-2 -1 D 1 2

Figure 4.48: Conducting strip image constructed using a bandwidth = 22.13%.

102

 

-2 -1 U 1 2

Figure 4.49: Conducting strip image constructed using a bandwidth ; 15.76%.

103

 

-2 -1 0 1 2

Figure 4.50: Conducting strip image constructed using a bandwidth = 12.26%.

104

 

-2 -1 U 1 2

Figure 4.51: Conducting strip image constructed using a bandwidth = 10.03%.

 

-2 -1 U 1 2

Figure 4.52: Conducting strip image constructed using a bandwidth = 8.49%.

106

 

-2 -1 0 1 2

Figure 4.53: Conducting strip image constructed using a bandwidth = 7.35%.

107

 

-2 -1 U 1 2

Figure 4.54: Conducting strip image constructed using a bandwidth = 6.49%.

108

 

-2 -1 D 1 2

Figure 4.55: Conducting strip image constructed using a bandwidth 2 5.81%.

4.2.2 Sampling Rate

In a frequency domain system, inadequate sample rate is the main cause of aliasing.
When the data is sampled at too large a frequency step size, the transformed time
domain data wraps around such that later events overlap with earlier events. This
effect produces an image with scattering centers that should be separated in space, but
instead overlap. To investigate the affect of aliasing on image quality, the frequency
step size of the simulated data is increased and the images are formed. The start and
stop frequency, 20 MHz to 16 GHz, remains the same in every case, but the number
of data points within this band is changed. The resulting images are shown in Figures

4.57-4.62, and the computed image radius is shown in Figure 4.56. As the frequency

109

step size increases, the image radius widens and the image becomes blurrier. At a step
size of 160 MHz (100 samples) the conducting strip is lost completely and therefore

no image radius is available.

 

25

IIIIIIIYU

N
I'I

Image Radius (cm)
2';
I .
Cl

05

 

 

 

IIUUIU'IIIIIUIII

l l I I I l 1 I I l l l
60 80 100 120

Sampling Rate (MHz)

-
1.
—
t
r

 

N
O
.b
0

Figure 4.56: Effect of sampling rate on image radius.

As seen in the figures, artifacts become more prominent as the sample step size
is increased. In Figure 4.58 the ghost of the strip is brighter, which could result in
false identification of another conducting strip. Figure 4.60 shows multiple artifacts
around the conducting strip creating background noise and producing a larger image
radius. The conducting strip is completely blends into the background in Figures 4.61

and 4.62.

110

0.5

-0.5

-1

 

~25

-2 -1 0 1 2

Figure 4.57: Conducting strip image constructed using a sampling rate of A f = 20

MHz.

111

0.5

-0.5

 

-2

-2 -1 0 1 2

Figure 4.58: Conducting strip image constructed using a sampling rate of A f = 40

MHz.

112

0.5

-1

 

-2

-2 -1 0 1 2

Figure 4.59: Conducting strip image constructed using a sampling rate of A f = 60
MHz.

113

     

-2 -1 0 1 2

Figure 4.60: Conducting strip image constructed using a sampling rate of A f = 80
MHz.

114

 

-2 -1 0 1 2

Figure 4.61: Conducting strip image constructed using a sampling rate of A f = 100
MHz.

115

 

-2 -1 0 1 2

Figure 4.62: Conducting strip image constructed using a sampling rate of A f = 120
MHz.

4.2.3 Signal-to—Noise Ratio

The same SNR issues seen with time domain data can occur in the frequency domain.
To investigate the effect of SNR on images created using frequency-domain data,
white Gaussian noise was added directly in the frequency domain using the program
WaveCalc. In this case, the SNR is defined as the ratio of signal power to noise power
(as opposed to a ratio of amplitudes as in the time domain). Figure 4.63 shows the

dependence of image radius on SNR.

116

 

S
I
III

d
0
III

0)
I l I I 1

Image Radius (cm)
0)
. ., . .

D

.. ..1...i...1...|...1...i...|.m.1
44 42 4O 38 36 34 32 30 28 26

SNR (dB)

I r T T '

 

 

 

Figure 4.63: Effect of the SNR on image radius.

This graph shows that image radius increases with decreasing SN R. Some discrep-
ancies exist wherein the image radius decreases with the decrease of SNR. This may
result from a wrong threshold value or a rogue peak of intensity showing up in the
noise. Figures 4.64 - 4.73 show the images created with various values of SNR. As

expected, the presence of noise makes it difficult to locate the strip.

117

 

-2 -1 0 1 2

Figure 4.64: Conducting strip image constructed using a SNR = 44 dB.

118

 

-2 -1 U 1 2

Figure 4.65: Conducting strip image constructed using a SNR = 42 dB.

119

0.5

-1.5

 

-2 -1 0 1 2

Figure 4.66: Conducting strip image constructed using a SNR = 40 dB.

120

0.5

-1.'5

 

-2 -1 U 1 2

Figure 4.67: Conducting strip image constructed using a SNR = 38 dB.

121

 

-2 -1 0 1 2

Figure 4.68: Conducting strip image constructed using a SNR = 36 dB.

122

0.5

-1.5

 

-2 -1 0 1 2

Figure 4.69: Conducting strip image constructed using a SNR = 34 dB.

123

 

-2 -1 0 1 2

Figure 4.70: Conducting strip image constructed using a SN R = 32 dB.

124

 

-2 -1 0 1 2

Figure 4.71: Conducting strip image constructed using a SNR = 30 dB.

125

0.5

-1.5

 

-2 -1 U 1 2

Figure 4.72: Conducting strip image constructed using a SN R = 28 dB.

126

0.5

 

-1.5

-2 -1 0 1 2

Figure 4.73: Conducting strip image constructed using a SNR = 26 dB.

4.2.4 Digitization

Digitization in a frequency domain system takes place directly in the frequency do-
main. To construct the images, the data is transformed into the time domain. An
important consideration is that in a frequency domain system, the reflection from the
target cannot be separated from the flash in the front of the barrier, as it can in a
time domain system (by using a time window). Thus, the target signal may have a
larger discretization error than with the time domain system.

Digitizing is done exactly as in the time domain, by assigning data to discrete
bins. The signals are moved into bins depending on the bin size described in (4.2),

and digitized based on (4.1). Figure 4.74 shows the image radius as a function of the

127

number of bits used to represent the data. Figures 4.75 - 4.85 show the images for

various numbers of bits used.

 

Image Radius (cm)
(0 h 01 O) \I on
IIUIUIIV'ITTT'IIII'ITTIIIIIT'T
E1

to
I
C]
D
C]

 

 

 

Bits

Figure 4.74: Effect of digitization on image radius.

A different type of blurring is seen than was observed with the digitization of the
time-domain data, and the target can just barely be distinguished with a discretiza-

tion of 4 bits.

128

0.5

 

-1.5

-2 -1 U 1 2

Figure 4.75: Conducting strip image constructed using the number of bits = 12.

129

0.5

 

-1.5

-2 -1 0 1 2

Figure 4.76: Conducting strip image constructed using the number of bits = 11.

130

0.5

 

-1.5

-2 -1 0 1 2

Figure 4.77: Conducting strip image constructed using the number of bits = 10.

131

0.5

 

-1.5

-2 -1 O 1 2

Figure 4.78: Conducting strip image constructed using the number of bits = 9.

132

0.5

 

-1.5

—2 -1 D 1 2

Figure 4.79: Conducting strip image constructed using the number of bits = 8.

133

0.5

 

-1.5

-2 -1 0 1 2

Figure 4.80: Conducting strip image constructed using the number of bits = 7.

134

0.5

 

 

-1.5

-2 -1 O 1 2

Figure 4.81: Conducting strip image constructed using the number of bits = 6.

135

0.5

 

-1.5

-2 -1 0 1 2

Figure 4.82: Conducting strip image constructed using the number of bits = 5.

136

0.5

-1

 

-1.5

-2 -1 0 1 2

Figure 4.83: Conducting strip image constructed using the number of bits = 4.

137

0.5

 

-1.5

   

-2 -1 0 1 2

Figure 4.84: Conducting strip image constructed using the number of bits = 3.

138

0.5

-1

 

-1.5

-2 -1 0 1 2

Figure 4.85: Conducting strip image constructed using the number of Bits = 2.

4.3 Two Conducting Strips

The parameter analysis is continued to compare the performance of the ultra—wide
band short—pulse and stepped frequency radar systems when imaging multiple targets.
Simulated data is computed at 21 different transmitter/ receiver positions, images are
created for different parameter sets and the sharpness of each image is analyzed. As

seen in Figure 4.86, the multiple targets used are two conducting strips.

139

 

l 21 different transmitter/receiver positions _
Equally spaced over 5 m _—‘l Z‘O-5m

...XX XX XX XX XX XX XX XX XX XX XX XX XX XX XX...

 

 

 

£0

 

 

z z=o.1524m
Bsr=6 0=.OOlS/m 8
£0 5 2:0
— -— =-O.8m

 

'Lo.2m* 0.4m io.2m*

Figure 4.86: Parameters of two conducting strip study

The program used to compute the field scattered by a single conducting strip
needs to be modified only slightly to accomodate scattering by the two conducting
strips shown in Figure 4.86. To compute the field scattered by the two strips, a largef
conducting strip of 0.80 m is first used to fill the Amn and bmn matrices of equations
(2.159) and (2.149). Then the entries corresponding to the 0.4 m gap between the
two strips are removed from these two matrices by replacing the stored values with
zeros. This models two strips of 0.2 m' in length, separated by a gap of 0.4 In.

The scattered field from the two conducting strips is found by first solving the
matrix equation (2.148). The solution is the complex current at each partition on
the two conducting strips. The total field is then determined at the receiver positions
by adding the scattered field from the known current on the conducting strips to the
incident field.

When the pulse width is altered in both the frequency and time domain, the
image radius changes for both conducting strips allowing for optimization of these
parameters in an actual prototype. Only the pulse width in the time domain and the

bandwidth in the frequency domain is examined in this section for multiple targets.

140

4.3.1 Pulse Width

With the geometry parameters of Figure 4.86, values of the computed image radii are

shown in Figure 4.87.

 

 

 

 

12' I l I I I
y.
10 "
E
3
- I
.3 8
'0 t
as
m .-
m 6'
m u-
m
g .
4 -
. I Left Strip
. ' I Right smp
+ '3
2 I
l J; #4 l l l I l l J L k l l J l I L L l l l L L I J mlm l l
0 0.5 1 1.5 2 2.5 3
Pulse Width (ns)

Figure 4.87: Effect of pulse width on image radius of two conducting strips.

This figure demonstrates that increasing the pulse width enlarges the image radii
which corresponds to an increase in the blurring of the strips. The left and right
strips both increase with similar tendencies as the pulse width expands. The images
of the different pulse widths are displayed in Figures 488- 4.97. When the pulse width
reaches 0.5 ns, shown in Figure 4.92, the conducting strips blend together becoming
one large conducting strip. Thus, a pulse width of less then 0.5 ns is needed resolve
these multiple targets. In Figure 4.94, where the pulse width is 1.5 ns, the image of
the conducting strips start to blur into the background. The size, position, intensity,
and shape of the strips cannot be determined when the pulse width is larger then 1.5

US.

141

[1.5

-1.5

 

-2

 

-2. 5

-2 -1 0 1 2

Figure 4.88: Two conducting strip image constructed using a pulse width of T = 0.05

HS

142

 

-2 -1 0 1 2

Figure 4.89: Two conducting strip image constructed using a pulse width of 7' = 0.10

US

143

0.5

   

-2.5

-2 -1 U 1 2

Figure 4.90: Two conducting strip image constructed using a pulse width of 7' = 0.15

ns

144

 

 

-2 -1 U 1 2

Figure 4.91: Two conducting strip image constructed using a pulse width of 7' = 0.25

US

145

0.5

 

-2.5

-2 -1 U 1 2

Figure 4.92: Two conducting strip image constructed using a pulse width of 7' = 0.50

US

146

0.5

 

-2.5

-2 -1 0 1 2

Figure 4.93: Two conducting strip image constructed using a pulse width of T = 1.00

US

147

 

-2 -1 0 1 2

Figure 4.94: Two conducting strip image constructed using a pulse width of 7' = 1.50

US

148

0.5

-U.5

 

-2.5

-2 -1 0 1 2

Figure 4.95: Two conducting strip image constructed using a pulse width of 7- = 2.00

US

149

0.5

 

-2.5

-2 -1 U 1 2

Figure 4.96: Two conducting strip image constructed using a pulse width of T = 2.5

D8

150

 

-2 -1 0 1 2

Figure 4.97: Two conducting strip image constructed using a pulse width of T = 3.00

11S

4.3.2 Bandwidth

As stated before, increasing the bandwidth of a stepped-frequency radar system in—
creases the resolution of the generated images. The dependence of image radii 0n
bandwidth for the two conducting strips is shown in Figure 4.98. Since the con-
ducting strips are positioned symmetrically about the z—axiz, the image radii seen
in Figure 4.98 increase almost at the same rate. As the bandwidth gets smaller the

image radius increases causing blurriness for the targets.

151

 

a)
r fr'

I Left Strip
0 Right Strip

Image Radius (cm)
#- U'I O) \i
l l I I I l l l I I l I l I I—l—I l—I—r
-
I . - .

CD
I
.

 

 

 

_
2 n I I I L I l 1 n I n I I I I n l I

i 100 80 60 4o 20 I
Bandwidth(%)

Figure 4.98: Effect of bandwidth on image radius of two conducting strips.

Figures 4.99 - 4.108 show the decrease in image quality as the bandwidth increases.
Following the progression though each of the images, the conducting strips and the
barrier become progressively unclear. When the bandwidth is at 22.13%, as seen in
Figure 4.101, the conducting strips are nearly indiscernible from one another; this
makes locating the separate objects almost impossible. At a bandwidth of 8.81%,
as shown in Figure 4.108, the two conducting strips have almost completely blended

into the background.

152

0.5

 

-2.5

-2 -1 0 1 2

Figure 4.99: Two conducting strip image constructed using a bandwidth of 110.38%

153

0.5

 

-2.5

-2 -1 0 1 2

Figure 4.100: Two conducting strip image constructed using a bandwidth of 36.75%

154

   

-2 -1 0 1 2

Figure 4.101: Two conducting strip image constructed using a bandwidth of 22.13%

155

0.5

 

~25

-2 -1 0 1 2

Figure 4.102: Two conducting strip image constructed using a bandwidth of 15.76%

156

0.5

 

-2.5

-2 -1 0 1 2

Figure 4.103: Two conducting strip image constructed using a bandwidth of 12.26%

157

0.5

 

-2.5

-2 -1 0 1 2

Figure 4.104: Two conducting strip image constructed using a bandwidth of 10.03%

158

 

 

 

0.5

   

-2.5

-2 -1 0 1 2

Figure 4.105: Two conducting strip image constructed using a bandwidth of 8.49%

159

0.5

 

-2.5

-2 -1 0 1 2

Figure 4.106: Two conducting strip image constructed using a bandwidth of 7.38%

160

0.5

 

-2.5

-2 -1 0 1 2

Figure 4.107: Two conducting strip image constructed using a bandwidth of 6.49%

161

 

 

0.5

 

-2.5

-2 -1 0 1 2

Figure 4.108: Two conducting strip image constructed using a bandwidth of 5.81%

162

 

 

Chapter 5

Experiment

A time—domain laboratory system [16] was constructed to investigate whether the
image techniques used with simulated data in the parameter study can be replicated
in practice. The laboratory system is not a dedicated radar, but rather uses standard
instrumentation .

This time-domain data was acquired using the MSU reflectivity arch range. This
arch range consists of two horn antennas, which can be arbitrarily located on an arch
of 20 ft in diameter. As seen in Figure 5.1, the transmitting and receiving wideband
TEM-horn antennas remain stationary while the barrier and strip target are moved,
providing equivalent results to moving the transmit/receive antenna pairs, as done
with the simulated data, and thereby simulating a rail SAR. The transmitting and
receiving antennas are American Electronic Laboratories H-1498 TEM horns each

with a bandwidth of 2 to 18 GHz.

163

 

~—'.——q

 

 

 

 

Lens—s % Lens
Receiving Horn *y <—

Transmitting Horn

Figure 5.1: Geometry of arch-range scattering system at MSU.

5.1 Experiment Set up

A 4 ft by 8 ft barrier was constructed using 1 / 2 inch construction drywall mounted
on pine wood of dimensions 2 inches by 4 inches. The barrier is able to roll by using
small caster wheels. A conducting strip was constructed using a 2 in by 8 ft segment
of aluminum tape applied to a low-density styrofoam column supported 2 ft behind
the barrier. The time domain response is measured using a Hewlett Packard digital
sampling oscilloscope (DSO) model HP54750A. As seen in Figure 5.2, this device has a

time-domain—reflectometer/time-domain-transmission (TDR/TDT) plug-in module,

164

model HP54753A, which provides 20 GHz and 18 GHz channels for the oscilloscope.

Bounda

   

Transmitting
Horn

   
  

Receiving
Horn

  

Whom

 

Digital Sampllng

  

Pulse Gem

Figure 5.2: Time domain scattering measurement configuration.

This TDR/TDT device has an integrated step generator, which sends steps from
channel 3 with a rise time of 45 ps and an amplitude of 200 mV. This unit triggers
a step generator made by Picosecond Pulse Labs (PSPL), part number 4015B. A
remote pulse head (PSPL 4015RPH) follows the pulse generator and creates another
step. The step signal is sent into a PSPL 5208—DC pulse-generating network where a
pseudo-Gaussian pulse signal is generated. This pulse, shown in Figure 5.3, is applied

to the transmitting antenna.

165

 

 

WT1T

-0.05
E
a) -0.1-
U, r-
9 l
6 h
> -0.15 -

 

 

 

l I I I l I I l g I I l I 4I I I. I I I l

0 0.2 L 0.4 0.6 0.8 1
Time (ns)

 

 

Figure 5.3: Pulse generated by the PSPL 5208—DC pulse generating network.

The frequency spectrum of the generated pulse can be determined by examining
its Fourier transform, which is shown in Figure 5.4. It is seen that the pulse has a

bandwidth of approximately 20 GHz.

166

 

 

 

 

 

1.—
$0.8”
'c .
3
'5- .
E0.6-
< .
13'
50'4
II
02*
l-
l...l...l...l...l..rl.i.l...I1..l...
0 2 4 6 8 10 12 14 16 18

Frequency (GHz)

 

Figure 5.4: Spectrum of the generated pulse.

Although the antennas are positioned as close to one another as possible to simu—
late a monostatic set up, there still exists a bistatic angle of 9b = 14.20. The position

of the antennas with respect to the target is shown in Figure 5.5.

167

 

Conducting

W W Strip

fl

 

  
 

Styrofoam
Column

 

 

 

 

 

Drywall Constructed
Dielectric Slab

Figure 5.5: Location of transmitting and receiving antennas.

Channel 4 on the HP54753A is used to receive the signal by sampling the incoming
voltage. A time window of 10 ns at 1 ns / div is used with 1024 sample points, giving
a sample period of about 0.01 ns or an equivalent sample rate of 100 Gs/s. Note
that this is significantly faster than the minimum required sample rate of 5 Gs/s
identified in the parameter exploration section of this report. The sampled signal is
average within the D80 256 times to increase the SN R. The received signal is stored
in a memory location in the D80. Clutter is eliminated by also storing a measured
background signal with no barrier or target present, which is subsequently subtracted
from the target signal.

Target data was measured for 21 equally spaced barrier positions spread over a

distance of 1.77 m, as shown in Figure 5.6.

168

 

I 21 different transmitter/receiver positions
I Equally spaced over 1.77 m

 

 

_—_1 z=2.4m

 

 

...XX XX XX XXXX XX XX XX XX XXXX XX XX XX XX...

 

 

 

 

 

 

 

80 Zl z=10.1cm
Drywall
F =
80 y 0
'r 1 .21 m 1'
Z=-600m

l——5cm ———l

Figure 5.6: Experimental setup.

It is important to note that thermal drift of the instruments can disturb the

results. Therefore, the equipment should warm up for a minimum of one hour before

measurements are taken. Also, the equipment is vulnerable to electrostatic discharge.

To eliminate this concern, one should always be grounded while taking measurements.

5.2 Experiment Calibration

The scattering measurements stored in the D80 are affected by the transfer functions

of the antennas (HT(f), HR(f)), the lenses (HTL(f), HRL(f)), mutual interactions

between the antennas (H A( f )), interactions between the target and its surroundings

(H SC( f)), and arch-range clutter (HC(f)). Therefore, a calibration procedure is

necessary to recover the target response. A block diagram of the measurement system

[17] is shown in Figure 5.7.

169

 

 

 

Transmitting Transmitting

Source: Pulse
Antenna Antenna Lens

 

21

 

 

 

 

 

gilt; @Z] H50) Hscm EB

 

 

 

Scatter Mutual Arch-range
Interactions Clutter
Receiving
+ HRL(f) Antenna Lens

Receiving

Antenna HR(f)

R....... a w ......

Figure 5.7: Block diagram model of measurement system.

From Figure 5.7 the waveform of an unknown target is formulated as

am = E(f)HT(f)HR(f){HA(f)

+HTL(f)HRL(f) [Han + H500) + Ham] } + N(f) (5.1)

where Hg( f ) is the system transfer function being solved for and N (f) is the envi-
ronment noise. A second measurement is done to find the background response (with
the target absent) which is used to obtain Hg( f) The background measurement is

subtracted from the unknown target response. When the time between these mea-

170

surements is long, the subtraction process will worsen from small changes in the test

environment. The formulated response for the background is

Rb(f) = E(f)HT(f)HR(f) {HA(f) + HTL(f)HRL(f)HC(f)} + N(f)- (5-2)

Subtracting the background response from the target response thus gives

Raw) = S(f) [Han + HSCm] . (5.3)

where S ( f) is the system transfer function given by

3(f) = E(f)HT(f)HR(f)HTL(f)HRL(f). (5-4)

With this system transfer function the unknown target response can be found. A
theoretically known response is measured to determine S (f ) A conducting plat is
chosen as the known target (calibrator) because its theoretical response is —1 [9].

The background-subtracted plate response is given as

RC—bU) = 80) [Hg (f) + Hscml (5.5)

where Hg( f ) = —1 is the theoretical response of the plate. Absorbers have been
placed around the targets in the reflectivity arch range to reduce the mutual in-
teractions. Therefore, the mutual interaction term H SC( f ) can be assumed to be
negligible or can be eliminated using time gating. The system transfer function is

then found to be
_ RC_b(f)

Sf— , 5.6
H Hg“) ( >

171

and the unknown target response is found using

Hfgm = W. (57)

Substituting (5.6) into (5.7) results in

, _R._b(f>Hg<‘(f)
Hsm— RC—b(f)

 

(5.8)

To demonstrate the calibration process, the time-domain response of a background

subtracted conducting plate is shown in Figure 5.8.

 

.0

o

N
I

 

 

l

l’ ‘I I r

Measured Signal (V)
_o o
52 8
T

-0.06 -

 

 

 

 

'0-08'1.L.L...m..rjnn.4...1
2 4

Figure 5.8: Raw measured waveform for a conducting plate.

First, this signal is transformed into the frequency domain using a Fourier transform
and truncated between 2.0 — 18 GHz, which is the band specification of the anten-
nas. This spectrum is then divided by the theoretical plate response of —1, and the

resulting system transfer function is shown in Figure 5.9.

172

 

7x10 ‘2_

6x10"2_

'T'

5m”:-

4x10”:-

Amplitude

3x10'12 E
2x10'12 :-

1x10“:-

 

 

 

Time (ns)

Figure 5.9: System transfer function.

With the system transfer function obtained, the target response is found by apply-
ing (5.7). A gaussian modulated cosine is used to window the function with a width
of T = 012713 and centered at 10 GHz. Then, to complete the calibration process the

results are transformed back into the time domain.

5.3 Experiment Results

The image created using the experimental data is shown in Figure 5.10 . From this
figure it can be seen that the barrier has an approximate width of .1 m and the

conducting strip is located at about .6 m from the back of the barrier.

173

 

0.5

 

-2.5

-2 -1 0 1 2

Figure 5.10: Image of barrier and strip created using time—domain response at re-

sponses measured at 21 different transmitter/ receiver positions.

The image radius of this experimental image is 2.61 cm. Comparing this result
to those simulated with the FORTRAN program shows a promising concentration
of conducting strip intensities. This suggests that the time domain imaging system
can provide accurate images of a target located behind a lossy barrier, validating the

simulations results.

174

 

Chapter 6

Conclusion

In this thesis, a two-dimensonal canonical problem is established that allows the
affects of various radar parameters on an ultra-wide band short-pulse and a stepped
frequency radar system to be studied. The canonical problem consists of a line source
positioned above a lossy dielectric slab and a target located behind the slab. The
interrogation cylindrical waves produced by the line source impinge on the barrier
and interact with targets placed behind the barrier. The resulting scattered field is
analyzed in the frequency domain at specific discrete frequencies. The simulated data
acquired is similar to data resulting form a UWB stepped frequency radar system. A
time-domain system is established by using an inverse Fourier transform to convert
the data into the time domain.

Analyzing the performance of these two radar systems is achieved by comparing
images which are constructed using a simple scattering-center technique. Factors such
as pulse width, sampling rate, number of bits, and SNR are examined in the time
domain to determine their affect on the quality of the image. Bandwidth, sampling
rate, SNR, and number of bits are examined in the frequency domain. In addition, a
laboratory time-domain radar system was constructed using instrumentation to allow

a validation of the simulation results.

175

 

 

Tables 6.1 and 6.2 reveal the dependence of target image on various important
parameters. These results can be summarized by providing a set of minimal require-
ments for producing a good image, and these can be useful for designing an actual
radar system. It must be emphasized that these requirements are subjective, since
they depend on the desired image quality, and are relative to the barrier and target

parameters used in the simulation.

 

Pulse Width Sampling Rate SNR Digitization r

S 0.10 ns 2 5 Gs/s 2 15 dB 2 6 Bits

Table 6.1: Minimum system requirements for a time domain radar system.

 

Bandwidth Frequency Stepsize SNR Digitization

 

2 36.75 % S 40 MHz 2 40 dB _>_ 6 Bits

Table 6.2: Minimum system requirements for a frequency domain radar system.

In chapter 5, the theory was verified using experimental results. An actual tran-
sient pulse with a 2—18 GHz bandwidth was used for measurements to simulate a
realistic time domain radar system. From the results of the measurements, it is con-
firmed that given certain radar parameters, a target can be imaged behind a lossy

barrier.

6.1 TOpics for further study

Possible follow-studies include the analysis of different targets behind a dielectric
slab. An in-depth look at different shaped objects would be useful in designing a
radar system with optimized system parameters. In addition to examining different

objects, an increase in the number of targets located behind the barrier would provide

176

further insight into system performance. Next, a method for quantifying the signal-
to—clutter ratio and the ghosting which results from the multiple reflections from
the barrier and target could prove useful. Finally, with the time domain system it
would be helpful to look at the affects of time-jitter, and to alter the number of

transmitter/ receiver positions, on the image radius.

177

 

Appendix A: Fortran Code for
Scattered Electric Field

PRleRAM Electric-Field_Scattered
f***************************************************
! For Dielectric Slab

.1 Scattering by a strip using MoM— TM

.' Then Sending the wave back through slab *

f***************************************************
I

f DECLARING VARIABLES TYPES

INIPLICIT REAL=I<8 (A—H,O—Z)

PARAMETER (MATSIZE2512)
PARAMETER (LIMIT=4096 ,IENW=4*LIMIT ,KEY=4)
PARAMEIER (TEST=1.D—-13)
(I)1\/lPLEX*16 Anm(l\/IATSIZE,1\-1ATSIZE) , Bm(MATSIZE) , DIAG (MATSIZE)
OONIPLEXHe J, K, x
INTEGER IPVT(MATSIZE), IFLAG, NN, MM M, N, I, 11, AA
REAL1<8 GAIN/MA, PI, MU, EPSO, EPSR, EPS, FREQ, w, KO, v
REAhs OMEGA, A, B, RK, 2, H, T, SIG, zs, Yobs, Zobs, YS
REAus DELTA, YM, YN, ANSI, ANSQ, ANS3, ANS4, ANss, ANSG
REAL*8 ANS7, ANS8, Emag, Ephase, ANSTl, ANST2, ANST5

REAMS ANST6, ANST8, ANST7, WORKGENW), FINTR

178

REAL=I=8 FINTI, FINT, AINTR], AINTR2, AlNTIl, AINTI2

INTEGER IVVORK(LIMIT)

EXTERNAL FINTR, FINTI, FINT, AINTRl, AINTR2, AINTIl, AINTI2
(XXVIPLEXHG ZINT, Umn, M11111, Ben, d1, d2, CURR, mum

REAL WALL, OBS-Y, OBS_Z, SOURCE_Y, SOURCEZ, STRIP, WIDTH

REAL PARTIS, POINTS, ENDrREQ, FREQSTEP

100mm TO ALL or PROGRAM
(mam /NA1\.1E1/YI\'I, DELTA, M, N

(Imam /NAMI22/K0, K, OMHEA, CURR, MU, PI, J
(mix-m /NA1\-IEB/ 2, H, T, 28, YS, Y

(mam /NA1\1E}4/ ECURR, Yobs, Zobs

!DECLARING VARIABLES
PI=3.14159265358979323846D0
J = DCMPIX(O.D0,1.D0)

MU = 4.DO*PI*1.D—7

EPSO = 8.854D—12

W= .075D0

WIDTH 2W
T = .1524D0

WALL = T
ZS = —.8D0

179

STRIP = ZS
YS = 0.01D0

SOURCEY = YS
Zobs = .6524d0

OBS_Z = Zobs — T

Yobs = —0.00dO

OBS_Y = Yobs
EPSR = 6.D0
EPS = EPSO*EPSR
SIG = .OOIDO
FREQ = 20.d6

OPEN(10,FILE=’Electric_Field-Scattered_Job_11.DAT’,STATUS=’
UNKNOWN)

OPEN(14,FILE=’PARAI\*IEIERS.DAT’ ,STATUS=’UNKNOWN’)
OPEN(13 ,FILE=’Center-Current-big_strip .DAT’ ,STATUSz’UNKNOWN’)
partis = 20.
DO 1:1, 800
Print *, Freq
if (2*W*FREQ*20.d0/(3.D8) .LT. 10.(.10) then
NN=20
else
NN = 2*W*FREQ*20.d0/(3.D8)

end if

180

 

 

DELTA = 2 . D0*W/NN
z = —.8d0
H = .6524d0
SOURCEZ = H — T
Y8 = 0.01D0
OMEGA = 2.DO*PI*FREQ
K0 = OMEARSQRHEPSOHVIU)
X: DCMPIx(EPSR,—SIG /(OI\K3A*EPSO) )
K = K0*SQRT(X)
CURR = DCt-iPLX(1.D0,0.DO)

ESUM = 0 . D0

!DEFINING YM AND CREATING Em MATRIX
DO M=1, NN
YM = —w + (l\=I*1.DO—.5D0)*DELTA
Y = YM

CALL BENINT (FINTR,O.dO,10.d0*K0,ANSI
)

CALL BENINT (FINTI ,0. d0 , 10. dO*K0,ANS2
)

Bm(M) = —2.D0*ANSl — 2.DO*J*ANS2

ICREATING Amn MATRIX
Do N=M, NN
YN = —w + (N*1.D0—.5D0)*DELTA

181

."-f'.'_i Milk._-..Is.‘v I-

 

“4.3

 

ENDDO
ENDDO

DOM: 2, NN

CALL BENINT (AINTRI ,0 . D0 ,SQRT(K0) ,
ANS3)

CALL BENINT (AINT11,0.DO,SQRT(KO),
ANS4)

CALL BENINT (AINTR2,0.d0,10.

d0*k0 ,ANSS)

CALL BENINT (AINT12,0.D0,10.

d0*k0 ,ANSG)

E
u

2 . D0*ANS3+J *2.D0*ANS4
D2 = 2.D0*J*(ANSEH—J*ANS6')
an -—— D1 + D2

a = DELTA/2.D0

b = (M—N) *DELTA

rk = K0

call Intel (A,B,RK,ZINT)
Umn = —((MU*OMH}A)/PI)*ZINT

Amn(M,N) = Umn + an

DON=1,M—1

ENDDO

ENDDO

Amn(M,N) = Amn(N,M)

182

 

 

 

ISOLVING FOR THE CURRENTS

CALL IXDMPC (512 ,NN,Anm,IPVT,DIAG)

CALL SOLVEC (512 ,NN,Anm,IPVT,Bm)

lSending current produced waves back through slab

999

Y = Yobs

H = T — ZS

Z = T — Zobs
DOM=1, NN

YM = —w + (l\1*1.DO—.5DO)*DELTA
vs = 'Y1\1

CURR = DELTA*Bm(M)

CALL BENINT (FINTR,0.d0, 10.

d0*K0, ANS7)

CALL BENINT (FINTI,0.d0, 10.

d0*K0, ANSS)
EBUM -_- ESLM + 2.D0*(ANS7 + J*ANS8)
END DO
01 = FREQ/169
O2 2 DREAL(EBUM)
03 = DII\1AG(ESU1\1)
WRITE (10,999) 01,02,03
format (3(1x,e12.5))

PRINT *, CENTER CURRENT=, ABS(BI\I(I.d0+NN/2.d0))

183

 

 

write(13,999) 01, DREAL(Bl\1(1.d0-+NN/2.d0)), DIMAG(BI\1
(1.d0+NN/2.dO))

FREQ = FREQ + 20.D6
END DO
POINTS 2 I
ENDFREQ = FREQ
PREQSTEP = FREQ/POINTS

CLOSE (10)

close (13)

WRITE(14 ,*) ’SLAB THICKNESS = ’, WAIL

W'RITE(14,*) ’POISTION OF OBSERVATION POINT ALONG Y = ’, OBS_Y

WRITE( 14 ,*) ’DISTANCE OF OBSERVATION POINT FROM WALL = ’ ,
OBS-Z

WRITE(14 ,1:) ’POSITION OF SOURCE ALONG Y = ’, SOURCEY

WRITE(14 ,*) ’DISTANCE OF SOURCE POINT FROM WAIL = ’, SOURCEZ

WRITE(14 ,*) ’DISTANCE OF STRIP FROM WALL = ’, STRIP
WRHE(14,*)’W1DI'H OF STRIP = ’, WIDTH

WRITE(14 ,*) ’NUMBER OF PARTITION POINTS = ’, PARTIS
WRITE(14 ,*) NUMBER OF POINTS = ’, POINTS
WRITE(14,*) ’STOP FREQUENCY = ’, ENDFREQ

WRITE(14 ,*) ’STEP FREQUENCY = ’, FREQSTEP

close (14)

END

184

 

'.
‘ I
I‘D-awn- ._

 

111111111llllllllllllllllllllllIlllllllll111111111111111111111

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

REAL*8 FUNCTION Fintr (KY)
REAL*8 KY
EXTERNAL Fint
CI)MPLEX*16 Fint

Fintr = DREAL (Fint (KY))
RETURN

END

11111111lllIlIlllllllllllllllllllllllllllllllllllllll111111111

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

REAL>I<8 FUNCTION Finti (KY)
REAL>1<8 KY
EXTERNAL Fint
mMPLEXH6 Fint

Finti = DIMAG (Fint(KY))
RETURN

END

1111111l11111111111111111111111111111111!lllllllllllllllllllll

REAL>I<8 FUNCTION Aintr1(KY)
REAL*8 KY
EXTERNAL Aintl

(I)MPLEX*16 Aintl

185

 

1‘— s!“ .,

 

Aintrl = DREAL (Aint1(KY))
RETURN
END

11111111111111lllllllllllllllllll1111lllllllllllllllllllllllll

REAL=I=8 FUNCTION Ainti1(KY)
REAL>I<8 KY
EXTERNAL Aintl
(I)MPIEX*16 Aint1
Aintil = DIMAG (Aint1(KY))
RETURN

END

11111111111111lllllllllllllllllllllll1111111111111111111111111

REAL=I=8 FUNCTION Aintr2 (KY)
REAL*8 KY
EXTERNAL Aint2
(I)MPLEX*16 Aint2
Aintr2 = DREAL (Aint2(KY))
RETURN
END

111111111llll11111llllllllllllllllllllllllIIIlllllllllllllllll

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

REAL*8 FUNCTION Ainti2 (KY)

186

~‘. .‘sl' i-“ ‘35" LAY:

 

‘1" "an.

 

REAL>I<8 KY
EXTERNAL Aint2
CDMPLEXarl6 Aint2
Ainti2 = DIMAG (Aint2(KY))
RETURN

END

1111111Illlllll111lllllllllllllllllllllllll111111lllllllllllll

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

(X)MPLEX*16 FUNCTION Fint (KY)
CDMPLEXHG P, Q, J, K, CURR
integer M, N
REAL>I<8 MU, PI, Z, H, T, ZS, Y, OMEBA, KY, K0, YM, DELTA, YS
(IIVINUN /NA1\1E1/YM, DELTA, M, N
(ITWWIN /NAME2/K0, K, Ohm CURR, MU, PI, J
(DWV'KTN /NAME3/ Z, H, T, ZS, YS, Y
IF (KO*KO .LT. KY*KY) THEN
P = —1.D0*J*SQRT(KY*KY—K0*K0)
ELSE
P = SQRT(K0*K0—KY*KY)
END IF
Q = SQRT(K*K—KY*KY)
IF (DREAL (Q) .LT. 0.D0) THEN

Q = —1.DO*Q

187

 

 

 

 

ENDIF
Fint = —(J*OMEI}A*CURR*MU/(2.d0*PI))*(COS(KY*(

Y—YS) ) ) *EXP(J*P*(Z—H+T) ) /((—P*P/Q) *SIN (Q*T
)—Q* SIN (Q*T) +2.D0*J*P*COS(Q*T) )

END

1111Ill1111111111111111111111111111111llllllilllllllllllllllll

(DMPIEXH6 FUNCTION Aintl (TAU)
(I)MPLEX*16 P, Q, J, K, R, CURR
integer M, N

REAL*8 MU, PI, Z, H, T, ZS, OMEEA, KY, K0, TAU, DELTA, RAT,
YM, Y, YS

CDMMON /NAMEI/YM, DELTA, M, N
cmMiN /NAMEZ/K0, K, OMEGA CURR, MU, PI, J
COMMON /NAME3/ 2, H, T, 25, YS, Y
KY=KO—TAU*TAU
IF (K0*KO .LT. KY*KY) THEN
P = —J*SQRT(KY*KY—KO*KO)
ELSE
P = SQRT(K0*KO—KY*KY)
END IF
Q = SQRT(K*K—KY*KY)
IF (DREAL (Q) .LT. 0.DO) THEN
Q = —1.D0*Q

188

 

'.
1"” EH:— .

Ii“

 

ENDIF
IF (KY .EQ. 0.D0) THEN

RAT = 1.D0
ELSE

RAT = SIN(KY*DELTA/2.D0) /(KY*DELTA/2.D0)
END IF

R = (( -—P*P/Q) *SIN(Q*T)+Q*SIN (Q*T) ) /((—P*P/Q) *
SIN (Q*T)—Q*SIN (Q*T) +2.d0*J*P*COS(Q*T))

Aintl = —(MU>IOMEI}A*DELTA/(2.DO*PI) ) *(COS(KY*(
M—N) *DELTA) *RAT*R*EXP( J * 2.D0*P*ZS) ) /(SQRT(
KO+KY))
RETURN

END

lilliilllllllllllllllllllllllHllllIllIIIIIHIIIIIIIIHHIIIII

(DMPLEXHG FUNCTION Aint2 (TAU)
(X)MPIEX*16 P, Q, J, K, R, CURR
integer M, N

REAL=I<8 MU, PI, Z, H, T, ZS, OMEEA, KY, KO, TAU, DELTA, RAT,
YM, Y, YS

(13mm /NAI\/IE1/Y1\/1, DELTA, M, N
W /NAME2/K0, K, OMEGA, CURR, MU, PI, J
(IMMON /NAMEB/ 2, H, T, ZS, YS, Y
KY=TAU*TAUi—K0

IF (KO*KO .LT. KYakKY) THEN

189

 

 

 

P = -J*SQRT(KY*KY—KO*KO)
ELSE
P = SQRT(K0*KO—KY*KY)
END IF
Q = SQRT(K*K—KY*IQ’)
IF (DREAL (Q) .LT. O.D0) THEN
Q = —1.D0*Q
ENDIF
IF (KY .EQ. 0.DO) THEN
RAT = 1.D0
ELSE
RAT = SIN(KY*DELTA/2.DO) /(KY*DELTA/2.DO)
END IF

R = (Q*Q—P*P)/(—(P*P+Q*Q) —2.DO*P*Q*((1.DO+EXP
(-J *2.DO*Q*T) ) /(1.DO—EXP(-J *2.D0*Q*T) ) ))

Aint2 = —(I‘vIU*OMEI}A*DELTA/(2.DO*PI) ) *(COS(KY*(
M—N) *DELTA) *RAT*R*EXP(J *2.DO*P*ZS) ) /(SQ,RT(
KY+K0))
RETURN

END

HIHIIIHIIHIIHllllillllllillllllllllllllHlllllllIllllllIl

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

subroutine benint (Func ,mWN,UP, ans)

implicit real*8 (a—h,O—z)

190

 

 

 

parameter (limit=4096,1enw=4*limit ,key=5)
parameter(pi=3.14159265358979323846d0)
parameter (test=1.d—18)

real*8 work(1enw), K0

integer iwork(limit)

external Func

call HQAG (Func, down, up, 1.d—200, test, KEY, reSO,
ABSERR, NEVAL, IER, LIMIT, LENW, LAST, IWORK, WORK

)
ans 2 RESO
return

end

llllllllllllllllllllllllllIlllllllllllllllllllllllIlllllllllll

subroutine Intcl (a,bl ,rk , zint)

! computes integral from O to infinity of the function
! sin(ax)cos(bx)/(x*sqrt(rk‘2—x“2))

implicit real=k8 (a—h,o—z)

parameter (limit=4096,1enw=4*limit,key=4)
parameter(pi=3.14159265358979323846d0)
parameter (test=1.d—8)

rea1*8 work(lenw)

real*8 IA,IB,IC,Il,12

complex>k16 zint

191

 

,I
tD'f‘EYO ’_

 

 

integer iwork(1imit)
external Fint1,Fintg,Fintgb
common /IC1/ aa,bb

common /IC2/ beta

common /IC3/ rrk

! upfact determines upper limit of integration

! for 12

upfact = 6.d0

b = abs(b1)
beta = rk
aa = a

bb = b

rrk = rk

! compute IB

argl = (a+b)*beta

f1 = exp(—arg1)*Ei(arg1) + exp(arg1)*E1(arg1)
arg2 = abS(a—b)*beta

f2 = exp(—arg2)*Ei(arg2) + exp(arg2)*E1(arg2)
if ((a—b) .lt. 0.dO) f2 =—f2

IB 2 (1.dO/(4.d0*beta))*(f1+f2)

! compute IC

down = 0.d0

up = rk

192

mi. .

 

 

 

call [HAG(Fint1, down, up, 1.d—200, test, KEY,
NEVAL, IER, LIMIT, LENW, LAST, IWORK, W)

IC = resO

! compute IA
IA 2 IB—IC

! compute 11
down = 0.d0
up = sqrt(rk)

call UQAG (Fintg, down, up, 1.d—200, test , KEY,
NEVAL, IER, LIMIT, LENW, LAST, IWORK, \MBK)

II = 2.d0*res0
! compute 12
down = 0.d0
up = upfactakrk

call IXQAG (Fintgb, down, up, 1.d—-200, test, KEY,
, NEVAL, IER, LIMIT, LENW, LAST, IWORK, WUiK)

12 = 2.dO*reSO + IA
zint = dcmplx(11,12)
return

end

rea1*8 function Fint1(x)
implicit real*8 (a—h,o—z)

common /IC1/ a,b

193

resO, ABSERR,
resO, ABSERR,
resO, ABSERR

 

 

 

common /IC2/ beta
Fintl = sin(a*x)*cos(b*x)/(beta*beta+x*x)
return

end

real*8 function FWntg(t)
hnpflicit real*8 (a—h,o—z)
common /IC1/ a,b
common /IC3/ rk
x ==rk—t*t
if (x .ne. 0.dO) then
Fintg = Sin(a*x)*cos(b*x)/(x*sqrt(rk+x))
else
Fintg = a/sqrt(rk)
end if
return

end

real*8 function Fintgb(t)
hnpflicit rea1*8 (a—h,o—z)
common /IC1/ a,b

common /IC2/ beta

common /IC3/ rk

x = rk+t*t

194

 

 

 

Fintgb = sin(a*x)*cos(b*x)*(1.d0/x — (sqrt(x*x—rk*rk)/(x*x+
beta*beta)))/sqrt(x+rk)

return

end

Real*8 function E1 (X)

 

 

I
! Purpose: Compute exponential integral E1(x)
! Input : x — Argument of E1(x)

! Output: E1 — E1(x) (x > 0 )

!

 

 

IMPLICIT DOUBLE PRECISION (A—H,O—Z)
IF (X.EQ.0.0) THEN

E1=1.0D+300
ELSE IF (x.LE.1.0) THEN

E1=—DLOG(X) +((((1.07857D—3*X—9.76004D—3)*X+5.519968D
—2)*X—0.24991055D0)*X+O.99999193D0)*X—0.57721566D0

ELSE
ESl =(((X+8.5733287401D0) *X+18.059016973D0) *X
+8.6347608925D0) *X+0.2677737343D0
ES2=(((X+9.5733223454D0) *X+25.6329561486D0) *x
+21.0996530827D0) *X+3.9584969228D0
E1=DEXP(—X)/X*ESl/ES2
ENDIF
RETURN
END

real *8 function EI(X)

 

 

I
! Purpose: Compute exponential integral Ei(x)
! Input : x ———— Argument of Ei(x)

' Output: EI —— Ei(x) ( x > 0 )

I

 

 

IMPLICIT DOUBLE PRECISION (A—H,O—Z)
IF (X.EQ.0.0) THEN
EI=—1.0D+300'
ELSE IF (X.LE.40.0) THEN
EI=1.0DO
R=1.0D0
DO 15 K=1,100
R=R>i<K>kX/ (K+1.0D0) “=2
EI=EI+R

IF (DABS(R/EI).LE.1.0D——15) GO TO 20

15 CONTINUE
20 GA=O.5772156649015328D0
EI=GA+DLOG(X)+X*EI

ELSE
EI=1.0DO
R=1.0D0
Do 25 K=1,20
R:R*K/X
25 EI=EI+R

196

 

EI=DEXP(X) /x* E1
ENDIF
RETURN
END
!***********************************************************
SUBROU'TINE IXDMPC (nbig ,N,A,IPVT,DIAG) I5.
3 DMMPOSFS MATRIX A TO TRIANGULAR FORM
i ”H OOMPLEX VERSION HM
l
I

PROM: LINEAR ALGEBRA BY GILBERT SIRANG .I’

 

INTEGER N,IPVT(nbig)

INTEGER NMLI ,J ,K,KP1,M
(DMPLEXH6 A(nbig ,nbig) ,DIAG(nbig)
(I)MPLEX*16 P,T

!
IPVT(N) = 1

NMl = N—l

 

DO 60 K=1,Nl\11

KP1 = K+1
! FIND PIVOT P
M = K

Do 10 I=KP1,N
10 IF (ABS(A(I,K)) .GT. ABS(A(M,K))) M=I
IPVT(K) = M

197

IF (M NH K) IPVT(N) = —IPVT(N)
P = A(M,K)

A(M,K) = A(K,K)

A(K,K) = P

DIAG(K) = P

IF (P .EQ. 0.0) GOTO 60

! COMPUTE MULTIPLIERS

20 Do 30 I=KP1,N

30 A(I,K) = —A(I,K)/P

! MERCHANGE RQWS AND (DLUMNS
Do 50 J=KP1,N

T = A(M,J)

A(M,J) = A(K,J)

A(K,J) = T

IF (T .EQ. 0.0) GOTO 50

Do 40 I=KP1,N

A(I,J) = A(I,J) + A(I,K)*T

40 CONTINUE

50 CONTINUE

60 CONTINUE

70 DIAG(N) = A(N,N)*IPVT(N)

RETURN

198

 

A:

 

 

END
!***********************************************************
SUBROUTINE SOLVEC (nbig ,N,A,IPVT,B)

I SOLVES MATRIX HQN AX=B WITH A DEIDMHEED BY DCDMPC
INTEGER N,IPVT(nbig)

INTEGER NMl,K,KB,KP1,KI\*'11,M, I F231
(DMPIEXHO A(nbig , nbig) ,B(nbig) 3" I.

(DMPLEX* 16 S

 

I FORWARD ELIMINATION I '

IF (N EC. 1) GOTO 30

 

NMI = N—l
Do 10 K=1,NMl
KP1 = K+1

M = IPVT(K)
S = B(M)
B(M) = B(K)
B(K) = S

Do 10 I=KP1,N
10 8(1) = B(I) +A(I,K)*S
! BACK SUBSTITUTION

Do 20 KB=1,NMl

KMl =N—KB

K = KMl+1

199

 

w
A
FF!
A
v

H

B(K)/A(K,K)

DO 20 1:1,KMl

20 8(1) =B(I) +A(I,K)*S

30 B(l) =B(1)/A(1,1) a
RETURN
END

 

LI:

 

200

 

Appendix B: Fortran Code for
Incident Electric Field

PR(X}RAM E1ectric-Field_Incident
I***************************************************
! For Dielectric Slab

! Scattering by a strip using MOM— 'IM

! Then Sending the wave back through slab *
!***************************************************

!DECLARING VARIABLES TYPES

IMPLICIT REAL*8 (A—H,O—Z)

PARAMETER (MATSIZE2512)

PARAMETER (LIMIT=4096 ,IEVW=4*LIMIT ,KEY=4)

PARAMETER (TEST=1.D—13)

(DMPLEX* 16 Amn(IVIATSIZEJV’IATSIZE) , Bm(MATSIZE) , DIAG(MATSIZE)
CDMPLEXH6 J, K, X

INTEGER IPVT(MATSIZE), IFLAG, NN, MM, M, N, I, II, AA
REAL*8 GAMMA, PI, MU, EPSO, EPSR, EPS, FREQ, W, KO, OMEEA
REAL*8 A, B, RK, Z, H, T, SIG, ZS, Yobs, Zobs, YS, Y
REAL*8 DELTA, YM, YN, ANSI, ANS2, ANSB, ANS4, ANS5, ANSO

REAL>I<8 ANS7, ANS8, ANSQ, ANSlO, ANSll, ANSI2, Emag, Ephase

201

 

 

REAL=I=8 “(Ei((IENW) , IINTRl, IINTII , IINTl, IINTR2, IINTIZ,
IINT2

INTEGER IWORK( LIMIT)
EXTERNAL IINTRI, IINTII, IINTl, IINTR2, IINTI2, IINT2

(DMPLEXHG ZINT, Umn, an, Ben, d1, d2, CURR, IBUM, TOTAL,
TOTALl, TOTAL2

real*4 01,02,03

ICIMVON TO ALL OF PROGRAM
(mam /NAME1/YM, DELTA, M, N

COMMN /NAME2/K0, K, OMEGA, CURR, MU, PI, J
COMMON /NAMEB/ 2, H, T, ZS, YS, Y

(13mm /NAME4/ ECURR, Yobs, Zobs

IDECLARING VARIABLES
PI=3.14159265358979323846D0
J = DCMPLX(0.D0,1.DO)

MU = 4.D0*PI*1.D—7

EPSO = 8.854D—12

W: .075D0
ZS 2: ——.8DO
YS = 0.01D0

= .1524D0

Zobs = .6524d0
Yobs = —0.00d0

202

 

 

.d '
‘mnl .

 

EPSR = 6.D0
EPS = EPSO*EPSR
SIG = .OOIDO

FREQ: 20.d6

OPEN(10,FILEz’Electric-Field-Incident-Job_11.DAT’,STATUS=’
UNKNOWN)

DO 1:1, 800
Print *, Freq
if (2*VV*FREQ*40.dO/(3.D8) .LT. 10.d0) then
NN=20
else
NN = 2*W*FREQ*20.d0/(3.D8)
end if

DELTA = 2 . DOarW/NN

Z = Zobs
H = .6524d0
Y8 = 0.01D0
Y = Yobs

OMEGA = 2.D0*PI*FREQ
K0 = Oh'IEEA*SQRT(EPSO*I\/M)

X: DCMPLX(EPSR, —SIG / (OMEEMEPSO) )
K = K0>I=SQRT(X)

CURR = DCIVIPLX(1.D0,0.DO)

203

 

total = 0.D0
CALL BENINT (IINTR1,0.D0,SQRT(K0) ,ANSll)
CALL BENINT (IINTIl ,O.DO,SQRT(KO),AN812)
CALL BENINT (IINTR2,0.d0,10.DO*KO,ANSQ)
CALL BENINT (IINT12 ,0.D0,10.D0*K0,ANSIO)
TOTALI = 2.DO*(ANS11+J*AN812)
TOTAL2 = 2.DO*J*(ANS9+J*ANSIO)

TOTAL = TOTALl + TOTAL2

 

01 = FREQ/1.d9
02 = DREAL(TOTAL)
03 = DIMAG(TOTAL)

WRITE (10,999) 01,02,03
999 format (3(1x,e12.5))

FREQ = FREQ + 20.D6
ENDDO
CLOSE (10)
close (11)

END

llllllliillIHHIIHIIHIlIII!IllllllllllllllllllllllIllllllll

REAL*8 FUNCTION Iintrl (KY)
REAL*8 KY

EXTERNAL Iintl

204

(DMPIEX*16 Iintl

Iintrl = DREAL (Iint1(KY))
RETURN
END

llllHllllllllIHlllHllllllIlllllIIHIIIIIHIIIIIllllllllllll

REAL*8 FUNCTION Iintil (KY)
REAL*8 KY
EXTERNAL Iintl
(X)MPLEX*16 Iintl

Iintil = DIMAG (Iint1(KY))
RETURN
END

lllllllllllIllllllllllllllllIlllllllllllllllllIlllllllllllllll

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

REAL*8 FUNCTION Iintr2 (KY)
REAL>|<8 KY
EXTERNAL Iint2
(DMPLEXIc16 Iint2

Iintr2 = DREAL (Iint2(KY))
RETURN

END

llllllllllllllllllllllllllllIlllllllllllilllllllllllllllllllll

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

ly'l'icf '33".

"". H z

 

 

REAL>I=8 FUNCTION Iinti2 (KY)
REAL>I=8 KY
EXTERNAL Iint2
(DMPLEXHG Iint2

Iinti2 = DIMAG (Iint2(KY))
RETURN

END

IlllllllllllllllllllllllllllllIllIlIllllllllllllHlllllllHlll

(DIV‘IPLEXHG FUNCTION Iintl (TAU)
(I)MPIEX*16 P, Q, J, K, CURR, R
integer M, N

REAL*8 MU, PI, Z, H, T, ZS, Y, OMEGA, TAU, KY, K0, YM, DELTA,
YS

(II/mm /NAME1/YM, DELTA, M, N
(DMNU‘I /NAME2/KO, K, OMEGA, CURR, MU, PI, J
(DI/mm /NAI\IE3/ Z, H, T, ZS, YS, Y
KY=K0—TAU*TAU
IF (K0*K0 .LT. KY*KY) THEN
P = -1.DO*J*SQRT(KY*KY—KO*KO)
ELSE
P = SQRT(K0*K0—KY*KY)
END IF
Q = SQRT(K*K—KY*KY)

206

 

 

IF (DREAL (Q) .LT. 0.DO) THEN
Q = —1.D0*Q
ENDIF

R = (Q*Q—P*P)/(—(P*P+Q*Q) —2.DO*P*Q* ((1 .DO+EXP
(—J *2.D0*Q*T) ) /(1.D(}-EXP(-—J *2.DO*Q*T) ) ))

 

Iintl = —(CURR*MUICMFI;A/(2.D0*PI))*COS(KY*(Y— E

YS) ) *(EXP(—J*P*ABS(Z—H) )+R*EXP(—J*P*(Z+H €..
—2.DO*T)))/(SQRT(K0+KY)) '

RETURN i .

lllllllllllIllllllIllllIllIIlllllllllllllllllllllllillllllllll

(DMPIEXHG FUNCTION Iint2 (TAU)
CI)MPIEX*16 P, Q, J, K, CURR, R
integer M, N

REAL*8 MU, PI, Z, H, T, ZS, Y, OMEGA, TAU, KY, K0, YM, DELTA,
YS

(IJMMON /NAME1/YM, DELTA, M, N
COMMON /NAME2/KO, K, OMEGA, CURR, MU, PI, J

COMMON /NAME3/ Z, H, T, ZS, YS, Y

 

KY=TAU>I<TAU+K0

IF (K0*KO .LT. KY*KY) THEN

P = —1.DO*J*SQRT(KY*KY—KO*K0)
ELSE

P = SQRT(K0*KO—KY*KY)
END IF

Q = SQRT(K*K—KY*KY)

207

IF (DREAL (Q) .LT. 0.D0) THEN
Q = —1.DO*Q

EN DIF

R = (Q*Q_P*P)/(—(P*P+Q*Q) —2.D0*P*Q
*((1.D0+EXP(—J*2.D0*Q*T))/(1,D0_
EXP(-J*2.DO*Q*T) ) ))

Iint2 = —(CURR*l\/IU*OMEEA/(2.D0*PI))* E
COS(KY*(Y—YS) ) *(EXP(—J*P*ABS(Z—H))
+R=I<EXP(—J*P* (Z—I—H—2.D0*T) ) ) /(SQRT(
K0+KY))

 

END

lllllllllllIIllIIIIIIIIll!IIHIIIIIHIIHIIIIIIIllllllllllllll

subroutine benint (Func ,II)WN,UP, ans)

implicit real*8 (a—h,o—Z)
parameter (limit =4096,lenw=4*limit ,key=5)
parameter(pi=3.14159265358979323846d0)
parameter (test=1.d—15)
real*8 work(lenw), K0
integer iwork(limit)
external Func
call IXQAG (Func, down, up, 1.d—-200, test, KEY, resO,

ABSERR, NEVAL, IER, LIMIT, LENW, LAST, IWORK, WORK
)

ans = RESO

 

208

return

end

 

209

 

Appendix C: Matlab Imaging Code

clear all

WWMW/WVMVVMWWWWMfWJfoWW/WWWMI
%Loading Files into Matrix %
VWWWWKMVM R“ I? (57' 67' ((3/7 (5737" c‘fuMWW/VIW 9672“ r7“ (‘MWQVI 9‘ 67-" M" WWWIZI

%Tau from wavecalc

 

tau = 0.12;
%Engery
Eng = 1/tau;

load Electric-Field-Scattered_Job_01.dat;
load Electric-Field-Scattered_Job-02.dat;
load Electric-Field-Scattered_Job_03.dat;
load Electric_Field-Scattered-Job_O4.dat;
load Electric-Field_Scattered_Job_05.dat;
load Electric_Field_Scattered_Job_06.dat;
load E1ectric-Field_Scattered-Job_07.dat;
load Electric_Field_Scattered-Job_08.dat;
load Electric-Field_Scattered-Job_09.dat;
load Electric-Field_Scattered_Job-10.dat;

load Electric-Field-Scattered_Job_11.dat;

210

 

load Electric-Field-Scattered-Job_12.dat;
load Electric-Fie1d-Scattered_Job_13.dat;
load Electric_Field-Scattered-Job_14.dat;
load Electric-Field-Scattered_Job_15.dat;
load Electric-Field-Scattered-Job_16.dat;
load Electric-Field-Scattered_Job_17.dat;
load Electric-Fie1d-Scattered_Job_18.dat;
load Electric_Fie1d_Scattered_Job_19.dat;
load ElectriC_Field_Scattered_Job_20.dat;
load Electric-Field-Scattered-Job_21.dat;
load Electric-Field-Incident-Job_allT.dat;
ES(: ,1) = Electric_Field_Scattered_Job_01(:,2)/Eng;
ES(: ,2) = Electric-Field_Scattered-Job_02 (: ,2)/Eng;
ES(: ,3) = Electric-Field_Scattered-Job_03 (: ,2)/Eng;
ES(: ,4) = E1ectric-Field_Scattered_Job_04 (: ,2)/Eng;
ES(: ,5) 2 Electric-Field-Scattered-Job_05 (: ,2)/Eng;
ES(: ,6) = Electric-Field-Scattered_Job_06 (: ,2)/Eng;
ES(: ,7) = Electric-Field-Scattered-Job-07(2 ,2)/Eng;
ES(: ,8) = E1ectric-Field_Scattered_Job_08 (: ,2)/Eng;
ES(: ,9) = Electric-Field-Scattered_Job_09 (: ,2)/Eng;
ES(: ,10) = Electric-Field_Scattered_Job_10 (: ,2)/Eng;

ES(: ,11) = Electric-Field-Scattered_Job_11(:,2)/Eng;

211

 

 

ES(:
ES(:
ES(:
ES(:
ES(:
ES(:
ES(:
ES(:
ES(:
ES(:
El(:
E1(:
E1(:
El(:
E1(:
E1(:
EI(:
El(:
E1(:
E1(:
E1(:
E1(:

EI(:

,13)

Electric-Field-Scattered_Job_12 (:,
Electric_Fie1d_Scattered_Job_13(:,
E1ectric-Field-Scattered_Job_14 (:,
Electric-Field-Scattered-Job_15 (:,
Electric_Field_Scattered-Job_16 (:,
Electric-Field-Scattered_Job_17(:,

Electric-Field-Scattered_Job_18 (:,

2)/Eng;
2)/Eng;
2)/Eng;
2l/EHS;
2l/Eng;
2)/Eng;

2)/Eng;

Electric-Field-Scattered..Job-19 (: ,2) /Eng;

Electric-Field-Scattered-Job-20 (:,

Electric_Field_Scattered_Job_21(:,

2)/Eng;

2l/Eng;

Electric-Field-Incident_Job_allT (: ,2) /Eng;

Electric-Fie1d-Incident_Job_allT(: ,

E1ectric-Field-Incident_Job-allT (:,

Electric-Field-Incident_Job_allT (:,

Electric-Fie1d-Incident_Job_allT (:,
Electric-Field_Incident_Job-allT (:,

Electric_Field_Incident-Job_allT (:,

Electric_Field_Incident_Job_allT (:,

Electric_Field_Incident_Job_allT (:,

Electric-Field_Incident_Job_allT (:
Electric-Field-Incident-Job-allT (:
Electric_Field_Incident_Job_allT (:

Electric-Field_Incident-Job-allT (2

212

2”Eng;
2l/Eng;
2)/Eng;
2l/Eng;
2)/Eng;
2l/Eng;
2)/Eng;

2l/Eng;

,2)/Eng;
,2)/Eng;
,2)/Eng;

,2)/Eng;

 

 

EI(:,14) = Electric-Field-Incident_Job-allT(:,2)/Eng;
EI(:,15) = Electric_Field_Incident_Job_allT(z,2)/Eng;
EI(:,16) = Electric_Field_Incident_Job-allT(:,2)/Eng;
EI(:,17) = Electric-Field_Incident_Job_a11T(:,2)/Eng;
EI(:,18) = Electric-Field-Incident_Job_allT(z,2)/Eng;
EI(:,19) = Electric-Field-Incident-Job_allT(z,2)/Eng;
EI(:,20) = Electric-Field_Incident_Job_allT(z,2)/Eng;
EI(:,21) = Electric-Field_Incident_Job_allT(:,2)/Eng;
WWWWM/‘MVI‘fIVJOWo/V’WWWMWWWWWWWWWWWWWMWQM
(3W1 11' 1 t1 ‘al I Z 8 ' c'W%%fW%WflWWWWWoWWWflXWuMM/WWE
%speed of light

c=3E8;
%speed of wave in wall

v=c/sqrt(6);
%Time shift

shift 2 1E—9;
%End of file

last=32768;
%Time step from E—file

delta_t=1.52E—l2;
%Number of points at end of file to Shift

NS = shift/delta_t;

NS = round(Ns);

213

%position in the z—direction of the tramsitters
22.65;

%determines the gride Size
delta=.01;

%Front side of wall in Z-direction
wall-ZA=.1524;

%Back side of wall in z—direction
wall_ZB=0;

%size of display window in Z-direction
zp_low=—2;

Zp_high=.5;

%size of display window in y—direction
yp_low=—2.6;
yp-high=2.6;

%Number of gride Points in z and y direction
Zp_delta=((Zp_high-—zp_low)/delta)+1;
yp-delta=((yp_high—yp_low)/delta)+1;

%First position of transmitter from E—field file
y-start=—2.5;

%Incrementing Size of transmitter position from E—Field file
y_delta=.25;

%number of transmitter positions

214

 

 

num_files=21;

%position somewhere between wall and object
threshold-dist=—.5;

%Position of sources to front of wall for thresholding
threshold-dist2=.2;

%incrementing size of the threshold l
threshold-delta-z=((threshold_dist-Zp_low)/delta)+1;

%number of gride points for the threshold positions

 

threshold-delta_z2=((threshold-dist2—Zp_low)/delta)+1; E
%Fixing ES and Bi files to include shift 9
for m:1:1:num_files;

for i=1:1:Ns;

ESO(i ,m) = ES(last—Ns+i ,m);
end

for i=Ns+1:1:last;

ESO(i ,m)=ES(i—Ns,m);

end

for i=1:1:Ns;

EIO(i ,m) = EI(last—Ns+i ,m);
end

for i=Ns+1:1:last;

EIO(i ,m)=EI(i—Ns,m);
end

215

 

end
%initializing value
for n=1:1: Zp_delta*yp_delta;
value(n)=0;
end
fM/W/WFA/éMMWMMWWVMWVWVMJWMWVWWMWFIRM/WAX)
(figfiffifigéfzidIWWfl/MXWWMKXVMWMAWJI/MWEWVWWVWVJWXVJE
y=y—Start;
n=1;
for m=1:1: nuIn_files;
for Zp=Zp_low : delta : Zp_high;
for yp=yp_low : delta: yp_high ;
dist=sqrt ((y—yp) ”2+(Z—Zp) A2);
if (Zp < wall-ZB)
thetazasin ((abs(zp—z) )/dist );
dist-A=abs(wall_ZA—Z )/sin ( theta);
dist-B=abs ( wall-ZB—wall_ZA )/sin ( theta);
dist-C=abs (Zp—wall-ZB )/Sin ( theta);

time = (2*dist-A/c)+(2*dist_B/v)+(2*dist-C/c)+
shift;

elseif (zp < wall_ZA)
thetazasin((abS(Zp-Z))/dist);
dist-A=abs(wall-ZA—Z)/Sin ( theta);

216

 

 

 

dist-B=abS(Zp—wall_ZA)/sin(theta);
time : (2*dist_A/c)+(2*dist-B/v)+shift;
else
time=2*dist/c+shift;
end
p=tiIne/delta_t;
post=round(p)+2;

post-low=round (p—delta/c)+2;

 

post-high=round(p+delta/c)+2;
if (post-low >= last)
A=ESO( last ,m) ;
else
A=[ESO( post-high ,m) ,ESO(postJow ,m) ,ESO(post ,m) I;
end
E=max(abs(A) );
value (n)=value (n)+abs(E);
n=n+1;
end'
end
y=y+y_delta;
n=1;

end

217

 

R‘MWM/IfMMWflfiKfl/JflI/I‘I/(‘IVI‘I/(lWW/flf/MfflWflMflMflfWX/MXXflVWfi
%I n c i d e n t

%
‘V‘Mfl/I/‘o/‘I/‘W/I‘A/M/Effl/M/J/I‘WL/W‘Ifl/‘I/‘IM/‘I/MfI/M‘Wt/WIZ‘M/fl/MVMVEM/MK)
y: y - s t a r t ;
n = 1 ;

for m=1:1:num_files;
for Zp=zp_low:delta:zp-high;

for yp=yp-low:delta:yp-high;

dist=sqrt ((y—yp) 32+(Z—Zp) A2);

if (zp < wall-ZB)
theta=asin((abS(Zp—Z))/dist);
dist-A=abs(wall_ZA—z)/sin ( theta);
dist-B=abs(wall_ZB—wall_ZA)/sin ( theta);

dist_C=abS(zp—wall..ZB)/sin ( theta);

time = (2*dist-A/c)+(2*dist-B/v)+(2*dist_C/c)+
shift;

elseif (Zp < wall-ZA)
theta=asin((abS(Zp—z))/dist);
dist_A=abs(wall_ZA—z)/sin ( theta);
dist-B=abS(zp—wall_ZA)/sin ( theta);

time = (2*dist-A/c)+(2*dist_B/v)+shift;

else

time=2*dist /c+shift;

end

218

' ALT—In)";

 

 

p=time/delta_t;

post=round(p)+2;

post-low=round(p—delta/C)+2;

post-high=round(p—l—delta/C)+2;

if (post_low >2 last)
A=E10(last ,m);

else

A:[E10(post-high ,m) ,EIO(postJow ,m) ,EIO(post ,m) ];

end
E=max(abs (A) ) ;
value (n)=value (n)+abs (E);
n=n+1;
end

end

y=y+y_delta;

n=1;

end

%V%WVVWMXVMMKWWXWWVMWF%W%
%Quantification %
N bac k = 1 0;

fa c = 1 .5;

Int=0;

box-z-end=—.65;

219

box_y_start =—.23;
box_z_start =—.95;
box-y-end=.23;

box-area==(abs ( box-y-end~—box-y-start ) ) *(abs ( box-z_end—
box_Z_Start ));

grid_z_start=(abs(Zp_low—box_z_start)/delta)+1;
grid-z-end=(abs(Zp_low—box_z_end)/delta)+1;
grid-y-start=(abs(yp_low—box_y_start)/delta)+1;
grid-y-end=(abs(yp_low—box_y_end)/delta)+1;
n=zp_de1ta*(grid_z-start —1)+grid_y-start;
i=1;
for Zzgrid-z_start :1: grid_Z_end
for y=grid-y-start :1:grid-y-end
Box1( i )=value ((z—1)*yp_delta+y) ;
zi(i)=floor ((i—1)/(grid_y_end—grid_y_start+1))+1;
yi(i)=i-(Zi(i)—1)*(grid_y_end—grid_y_start+1);
i=i+1;
end
end
Box2=sort(Box1, ’descend ’) ;
iendzi—l;
sum_y=0;

sum_z=0;

220

sum_f=0;

sum_R=0;

sum-A=0;

for i=1:iend;
sum_y=sum_y+Box1( i )* yi ( i) ;
sum_z=sum_z+Box1( i )* zi ( i) ;
sum-f=surn-f+Box1( i) ;
y-0=sum_y/su1n_f;
Z_0=sum_Z/sum_f;

end

thresh=0.5*Box2(1);

for i=1:iend;

Di()i))=sqrt ((yi(i)—y-0)*(yi(i)—y-0)+(zi(i)—z_0)*(zi(i)—Z_0

Ai( i )=abS ((yi(i)—Y-0)*(Zi ( i l—Z-Oll;
if (Box1( i )>thresh)
sum_R=sum_R+Box1 ( i )*Di( i ) ;
sum-A=sum_A+Box1( i ) *Ai( i ) ;
end
end
RadiuS$um_R/sum_f;
Area=sum_A/sum_f;

v=sqrt (var (Boxl) );

221

 

fract1=0.99;
numl=0;
for i=1:iend
if (Box2(i) > fract1*mean(Box2))

num 1=numl + 1 ;

end
end
Box_roth=0;
for i =1:10

Box_roth=Box_roth+Box2( i );
end
Rothwell=Box-roth/1O
numl/iend
figure(2)
plot(Box2);
WWW/“WWWV5MMMMWMMWJWMWVMWMMMMWWMWWI/WWQ/“WM
$$$£fii§wmWMXWM/VM/MWWMZ‘W/M/é/MX/{MJ/‘I/WWMWEI
%for the stip
for m=1:1:threshold-delta_z*yp-delta;

threshold=0;

if (value(n) < threshold)

value(n)=0;
end

222

end
%for the slab

for n=round( threshold-delta_z*yp_delta) +1:1:round(
threshold-delta_z2*yp_delta);

threshold=0;
if (value(n) < threshold)
value(n)=0;
end
end
%for above Slab

for n=round(threshold-delta_z2*yp_delta) +1:1: yp_delta*
zp-de1ta

threshold=0;

if (value(n) < threshold)
value(n)=0;

end

end

ZWfiFiW/flWflWM/VMJMWflfi/wWWI/VMMMWWWWW/JE
for z=1:1:zp_delta;
for y=1:1:yp_delta;

F(Z ,y) = value(y+yp_delta*(z—1));

end

223

end

[y,z] = meshgrid(1:1:yp_delta ,1:1: Zp_delta);
figure(1)
pcolor(yp_low+delta*(y—1),Zp_low+delta*(z—1),F(: ,2) );
colorbar;

shading interp;

colormap hsv;

224

 

Bibliography

[1] Ferris Jr., D. D. and N. C. Curris, A survey of current technologies for through-
the-wall surveillance (TWS), proceeding SPIE, Vol. 3577, 6272, Jan. 1999.

[2] Ahmad, F. and M. Amin, A non coherent approach to through-the-wall radar
imaging, Proceedings of the 8th International Symposium on Signal Processing
and Its Applications, Vol. 2, 28 31, 2005.

[3] Taylor, J. D., Introduction to Ultra-wideband Radar Systems, CRC Press, 1995.

[4] Chen, F. C. and W. C. Chew, Time-domain ultra-wideband microwave imaging
radar system, Journal of Electromagnetic Waves and Applications, Vol. 17, No.
2, 313331, 2003.

[5] Maaref, N., Millot, P., and X. Ferri‘eres, Electromagnetic Imaging Method Based
on Time Reversal Processing Applied to Through-the— Wall Target Localization,
Progress In Electromagnetics Research M, Vol. 1, 5967, 2008.

[6] Soumekh, M., Synthetic Aperture Radar Signal Processing, John Wiley and sons,
Inc., New York, NY, 1999.

[7] Chan, Y. K. and V. C Koo, An introduction to synthetic aperture radar (SAR),
Progress In Electromagnetics Research B, Vol. 2, 2760, 2008.

[8] Charvat, G.L., A Low-Power Radar Imaging System, Ph.D. Thesis, Michigan
State University, 2007.

[9] D.K. Cheng, Field and Wave Electromagnetics Addison-Wesley Publishing Com-
pany, Inc, 1992.

[10] E.J. Rothwell and M.J. Cloud, Electromagnetics, Second Edition, CRC Press,
Boca Raton, Florida, 2009.

[11] R. Piessens, et. al., Quadpack: a subroutine package for automatic integration,
Springer-Verlag, Berlin, 1983.

[12] W. Press, B. F lannery, S. Teukolsky, and W. Vetterling, Numerical Recipes,
Cambridge University Press, Cambridge, 1986.

225

 

[13] N. Kinayman and MI. Aksun, Comparative study of acceleration techniques for
integrals andseries in electromagnetic problems, Antennas and Propagation Soci-
ety International Symposium, 1995.

[14] J Davis, Y Huang, SG Millard, and JH Bungey, Determination of dielectric
properties of insitu concrete at radar frequencies, International Symposium (NDT-
CE 2003), Non—Destructive Testing in Civil Engineering, 2003.

[15] http://netlib.org/quadpack/index.html.

[16] A. Kizilay, A perturbation method for transient multipath analysis of electromag-
netic scattering from targets above periodic surfaces, Ph.D. Thesis, Michigan State

University, 2000.

[17] J. E. Ross, III, Application of transient electromagnetic fields to radar target
discrimination, Ph.D. Thesis, Michigan State University, 1992.

226