Title: Array processing for composite wavefront decomposition
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00098196/00001
 Material Information
Title: Array processing for composite wavefront decomposition
Physical Description: xii, 233 leaves. : illus. ; 28 cm.
Language: English
Creator: Halpeny, Owen Simeon, 1940-
Publication Date: 1972
Copyright Date: 1972
 Subjects
Subject: Electric filters   ( lcsh )
Wave mechanics   ( lcsh )
Electric waves   ( lcsh )
Electromagnetic waves   ( lcsh )
Electric noise   ( lcsh )
Electrical Engineering thesis Ph. D
Dissertations, Academic -- Electrical Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
 Notes
Thesis: Thesis -- University of Florida.
Bibliography: Bibliography: leaves 229-232.
General Note: Typescript.
General Note: Vita.
 Record Information
Bibliographic ID: UF00098196
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000577349
oclc - 13958869
notis - ADA5044

Downloads

This item has the following downloads:

arrayprocessingf00halp ( PDF )


Full Text












Array Processing for Composite Wavefront Decomposition


By

Owen Simeon Halpeny












A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA IN PARTIAL
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY


UNIVERSITY OF FLORIDA
1972
































PLEASE NOTE:


Some pages may have

indistinct print.

Filmed as received.


University Microfilms, A Xerox Education Company































To Michael Owen Halpeny














ACKNOWLEDGEMENTS
The author wishes to express his sincere appreciation to Dr.

D. G. Childers, chairman of his supervisory committee, for his

invaluable and timely inspiration, suggestions and encouragement.

He also wishes to thank Dr. T. S. George and Dr. N. W. Perry, Jr.,

the other members of his supervisory committee for their guidance.

The valuable discussions with Miss Julia McCoy and other staff

members of the Visual Science Laboratory are gratefully acknowledged.

The author also wishes to especially thank Mrs. Carol Halpeny

who assisted in the coding of the programs, typed the manuscript,

prepared the drawings, and offered many valuable comments regarding

the organization of this dissertation.















TABLE OF CONTENTS


ACKNOWLEDGEMENTS...............................................

LIST OF TABLES.................................................

LIST OF FIGURES...............................................

ABSTRACT ..................................... ..................

CHAPTERS

I INTRODUCTION........................... ... .........

II MULTIDIMENSIONAL POWER SPECTRAL DENSITY ESTIMATION...

Introduction.......... ............................

One-dimensional Spectra.............................

Frequency-wavenumber Power Spectral Density
Estimation ................................. ....

Two-dimensional Spectra .............................

Three-dimensional Spectra............................

Implementation....................... .. ... ... ..

Results ................................... .....

Summary.........................................

III MULTIDIMENSIONAL DIGITAL FILTERING...................

One-dimensional Filtering...........................

Three-dimensional Filtering.........................

Implementation......................................

Angle Filtering..................................

Digital Angle Filter Implementation ................


Page

ii1

vi

vii

xi


1

6

6

7


19

22

29

37

41

64

65

65

79

85

86

87










Velocity Filtering.................................. 106

Digital Velocity Filter Implementation............... 107

Frequency Filter.................................... 112

Frequency Filter Implementation..................... 114

Summary........................................ 114

IV MULTIWAVE MAXIMUM-LIKELIHOOD ESTIMATION ............ 117

Single Wave Maximum-likelihood Estimation............ 117

Multiwave Maximum-likelihood Estimation.............. 129

Implementation ........... ........ ................... 133

Examples .................................... ....... 136

Summary ........................ ...................... 148
V COMPUTER PROCESSING OF SPATIO-TEMPORAL DATA......... 149

VI SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS............ 162

Multidimensional Power Spectral Density Estimation... 162

Multidimensional Digital Filtering.................. 163

Multiwave Maximum-likelihood Estimation.............. 164

Computer Processing of Spatio-temporal Data.......... 165

APPENDICES

A NONPARAMETRIC LEAST MEAN SQUARE ERROR ESTIMATION..... 168

B A BRIEF REVIEW OF MAXIMUM-LIKELIHOOD ESTIMATION...... 174

C DERIVATION OF THE TIME DOMAIN SOLUTION............... 178

D PROGRAM LISTINGS................................... 187

REFERENCES .............. ............................... ........ 229

BIOGRAPHICAL SKETCH........................................... 233















LIST OF TABLES


Table Page

1 Number of operations................................ 37

2 Listing of overlapping planewaves .................... 43

3 Grey levels.......................................... 43

4 Number of operations required for various data
array sizes................................. ...... 85

5 Approximate resolution for various bearings........... 95






























vi














LIST OF FIGURES


Figure Page

2.1 Fourier transform pair for a triangular window....... 13

2.2 Fourier transform pair for a rectangular window...... 13

2.3 Fourier transform pair for ideal sinusoid............ 14

2.4 Fourier transform pair for rectangular-truncated
sinusoid....................................... .. 14

2.5 Fourier transform pair for triangle-truncated
sinusoid...................................... 14

2.6 Relative execution times.......................... 18

2.7 Two-dimensional monochromatic spatio-temporal
function with frequency f0, wavenumber kx0, and
velocity vx... .................................... 24

2.8 Frequency-wavenumber spectrum of two-dimensional
monochromatic spatio-temporal function with
frequency fO, wavenumber kx0, and velocity vx ....... 25
2.9 Frequency-wavenumber spectrum for two-dimensional
wideband spatio-temporal function with
velocity vx0 ... ................. .... .. ......... 26
2.10 Frequency-wavenumber spectrum of three-dimensional
monochromatic spatio-temporal function with
frequency fo, wavenumber RO, and velocity V0......... 31

2.11 Frequency-wavenumber PSD estimation by Time Domain
Method....................................... ......... 33
2.12 Simplified flow chart for simulated data............. 38

2.13 FWNO output, f = 0.0 Hz............................. 44

2.14 FWNO output, f = 1.6 Hz............................. 45

2.15 FWNO output, f = 3.1 Hz............................. 46

2.16 FWNO output, f = 4.7 Hz............................. 47










2.17

2.18

2.19

2.20

2.21

2.22

2.23

2.24

2.25

2.26

2.27

2.28

2.29


FWNO

FWNO

FWNO

FWNO

FWNO

FWNO

FWNO

FWNO

FWNO

FWNO

FWNO

FWNO

FWNO


2.30 Example of positive and negative frequency terms.....

3.1 Minimum transition width low-pass filter
amplitude response.................................

3.2 Minimum transition width low-pass filter
impulse response................ ....................

3.3 Impulse response with zeros added ...................

3.4 Minimum transition width interpolated
amplitude response.................................
3.5 Hanning window transition low-pass filter
amplitude response.................................
3.6 Hanning window transition low-pass filter
impulse response ...................................
3.7 Hanning window transition interpolated
amplitude response ................................
3.8 Optimum 2-point transition low-pass filter
amplitude response ..................................
3.9 Optimum 2-point transition low-pass filter
impulse response ...................................

viii


output, f = 6.3 Hz.............................

output, f = 7.8 Hz............................

output, f = 9.4 Hz............................

output, f = 10.9 Hz.........................

output, f = 12.5 Hz............................

output, f = 14.1 Hz.............................

output, f = 15.6 Hz.........................

output, f = 17.2 Hz...................... ......

output, f = 18.8 Hz............................

output, f = 20.3 Hz...........................

output, f = 21.9 Hz...........................

output, f = 23.4 Hz............ ..........

output, f = 25.0 Hz............................










3.10 Optimum 2-point transition interpolated
amplitude response................................... 77

3.11 Generalized multichannel filter...................... 80

3.12 Generalized multidimensional filter................. 80

3.13 Frequency response for one wideband pass filter...... 84

3.14 Angle filter frequency-wavenumber response........... 88

3.15 Number array output of subroutine FA................ 91

3.16 Grey tone contour output of subroutine FA............ 92

3.17 Simplified flow chart for angle filtering exercise... 94

3.18 Estimated waveform using only one channel............. 96

3.19 Estimates of planewave number 1...................... 98

3.20 Estimates of planewave number 2..................... 99

3.21 Estimates of planewave number 3..................... 100

3.22 Estimates of planewave number 4..................... 101

3.23 Estimates of planewave number 5..................... 102

3.24 Estimates of planewave number 6..................... 103

3.25 Estimates of planewave number 7....................... 104

3.26 Velocity filter frequency-wavenumber response........ 108

3.27 Number array output of subroutine FV................. 110

3.28 Grey tone output of subroutine FV................... 111

3.29 Frequency filter frequency-wavenumber response....... 113

3.30 Frequency filter example, planewave number 1......... 115

4.1 General maximum-likelihood estimator, K=2............ 121

4.2 Maximum-likelihood estimator with zero cross-power... 123

4.3 Optimum maximum-likelihood estimator for noise
independence between sensors with common
spectral density.................................... 124










4.4 Maximum-likelihood estimator form suitable for
numerical calculations............................. 126

4.5 General multiwave maximum-likelihood estimator....... 131

4.6 Simplified flow chart for maximum-likelihood
operation................. ......... .. ............ 134
4.7 Estimation comparison, SIR=2....................... 138

4.8 Estimation comparison, SIR=1........................ 139

4.9 Estimation comparison, SIR=0.5...................... 140

4.10 Estimation comparison, SIR=0.1..................... 141

4.11 Estimation comparison, SNR=0.5...................... 143

4.12 Estimation comparison, SNR=1.0..................... 144

4.13 Estimation comparison, SNR=5.0....................... 145

4.14 Multiwave estimation, wave number 2............... 146

4.15 Multiwave estimation, wave number 4 ................. 147

5.1 Scalp electrode layout (rear view).................. 150

5.2 Flow chart for PRELIM ............................... 152

5.3 Aliasing chart, 125 Hz sampling..................... 152

5.4 VER, white stimulus ................................. 154

5.5 VER, red-green stimulus............................. 155

5.6 VER, noise......................................... 156

5.7 FWNO output, 1.6 to 7.8 Hz, before filtering......... 157

5.8 FWNO output, 1.6 to 7.8 Hz, after filtering.......... 158

5.9 Delay and sum output, white stimulus................. 159

5.10 Delay and sum output, red-green stimulus............. 160

A.1 General least mean square error estimator............ 170

B.1 General statistical model............................ 174












Abstract of Dissertation Presented to the
Graduate Council of the University of Florida in Partial Fulfillment
of the Requirements for the Degree of Doctor of Philosophy

ARRAY PROCESSING FOR COMPOSITE WAVEFRONT DECOMPOSITION

By

Owen Simeon Halpeny

August, 1972
Chairman: Dr. D. G. Childers
Major Department: Electrical Engineering

This study considers the problem of detecting and estimating

multiple planewaves in a noise environment. The data are assumed

to be available in array form, i.e., several spatially distinct time

records of the same process. The solution which is pursued is

computationally efficient, but generally suboptimal. However, under

some.conditions, the solution is optimum in the maximum-likelihood

sense.

The determination of the number of planewaves and their respective

vector velocities is accomplished by using the multidimensional

frequency-wavenumber power spectral density function. This technique

avoids the complexity of more elaborate schemes and yields usable

results for sufficiently high signal-to-noise ratios. The multiwave

problem is reduced to a sequence of single wave estimation problems

by using frequency-wavenumber digital filters which pass a single

planewave. The multidimensional fast Fourier transform is used for

both the spectral estimation and digital filtering.











After filtering, maximum-likelihood estimates are made for each

wave. The estimates are compared to the simpler delay and sum

estimates.

All of the techniques advanced in this study are tested using

simulated and actual data on a large scale digital computer. It is

found that the spectral density estimation and filtering can be

accomplished efficiently and effectively. By comparison, the general

maximum-likelihood estimators are computationally cumbersome and

offer little improvement over the delay and sum estimates.















CHAPTER I

INTRODUCTION


The analysis of spatio-temporal signals which are immersed in

noise is undertaken here because of the current interest in many fields

of both theoretical and practical considerations. Typical research

includes the radar and sonar problem of locating and classifying

targets that emit or reflect waves of energy. In radio astronomy the

purpose is to study electromagnetic radiation from different regions of

the universe. In seismology one is frequently interested in E1]

differentiating between earthquakes and nuclear explosions.

The original motivation for this study is based upon the conjecture

that these sophisticated techniques may be of value in the study of

electroencephalographs (EEG) or brain waves. Accordingly, this study

is limited to techniques which have definite practical merit. The

signal model is a sum of planewaves which overlap in space and time.

The noise consists of both ambient channel noise and interfering

planewaves. It is generally assumed that the noise statistics are

available, but the signals are completely unknown a priori. Thus,

while this study is guided by a particular application, the formulation

of the problem is sufficiently general to permit application to many

other areas.










The basic quantity being studied is a wave. In general terms a

wave is simply a phenomenon or disturbance which moves. The shock wave

caused by an explosion or the tidal wave caused by an earthquake are

obvious examples. A planewave is simply a wave which travels in a

straight line and does not exhibit variations in any plane perpendicular

to the path of propagation. In the examples cited the wave travels away

from the source of disturbance. In complex biological systems this

simple relationship may not exist.

Mathematically a planewave may be defined as


s(t d/v)


where t is time, d is distance as measured from the spatial origin

along the direction of propagation, v is the velocity of the wave, and

s(t) is the waveshape or functional form of the wave at the spatial

origin. Using the conventional Cartesian coordinate system, this wave

becomes


s(t .4r)


where = (r ,r ,r ) is the position vector at which s is evaluated

and S is the inverse velocity vector


A monochromatic planewave in complex exponential form is given by










-*
-12in(ft-fi-r) = e-i2n(ft-'*7)


where the vector wavenumber is t = fl. Any spatio-temporal function

which is transformable may be expressed as a sum of monochromatic

planewaves.


s(t,-) = if S(f,)e-i21(ft-' ) df dl


where s(t,r) and S(f,t) form a multidimensional Fourier transform pair.

This provides an additional incentive for studying planewaves.

To determine the waveshape of a single wave of known vector

velocity, in a noiseless environment, one simply observes the wave at

any point in space and applies the appropriate delay (or advance) to

produce s(t). If a single planewave is spatially sampled using a

discrete array of sensors and each sensor adds white noise, an estimate

of the wave is obtained by aligning the signals on each of the channels

and summing the results. This is the simplest form of array processing.

If, in addition to the channel noise, noise exists in the form of an

interfering planewave, more complicated processing must be employed.

This is a common problem encountered in radar. Allen [2] provides a

particularly readable introduction to this subject and lists some of

the many references. This estimation problem is usually solved by

adding amplitude weighting to each channel before the delay and sum

operations. Design engineering judgment is usually applied to obtain

acceptable estimates under a variety of conditions.










To obtain the optimum estimates, it is necessary to apply the

concepts of decision and estimation theory. For most cases of practical

interest, the array processor then requires a different linear filter

applied to each channel before the delay and sum operation. The

filters generally depend upon the array geometry and noise statistics.

The resulting estimates are referred to as maximum-likelihood estimates.

This problem was first solved by Kelly and Levin E13. A closed form

solution was not obtained to the general problem of estimating the

waveshape and vector velocity. An integral equation must be solved by

iteration. The parameter being varied is the vector velocity. If the

vector velocity is known, a closed form solution is possible. Capon

[3] has successfully implemented this solution and has applied the

results to seismological data. In a later paper Capon [4] also

discusses the wave detection problem.

The next step is to consider the multiwave maximum-likelihood

estimation problem. For the special case where the number of signal

and noise planewaves is known, the solution is given by Schweppe [5].

Again, a closed form solution is not obtained. This time it is

necessary to vary all of the vector velocities simultaneously in order

to solve an integral equation by iteration. When compared to the one

wave problem, this is referred to as a multidimensional search. If the

vector velocities are known, closed form solutions are possible. An

implementation is given by Kobayashi and Welch [63 for two planewaves

of known vector velocities. The multiwave detection problem, i.e.,

determining the number of waves, has not received much attention in

the literature.










The general problem considered in this dissertation is to

determine the number, vector velocity, and waveshape of overlapping

planewaves in the presence of additive noise. A general optimum

solution is not found. Instead, a heuristic solution is presented

along with a complete working implementation scheme for large scale

computers. For the case where the number of waves and the vector

velocities are known, the solution is optimum.

The detection of waves and the estimation of the vector velocities

is accomplished heuristically by using the frequency-wavenumber power

spectral density function. A formal solution of this part of the

problem would require maximum-likelihood detection and an exhaustive

search algorithm. The general subject of frequency-wavenumber power

spectral density estimation is covered in Chapter II.

To reduce the multiwave problem to a succession of single wave

estimation problems, the frequency-wavenumber filters are employed,

which pass only a single planewave. This avoids the difficult multi-

dimensional search required by Schweppe's solution. The general

subject of frequency-wavenumber filtering is covered in Chapter III.

The last step is to form a series of single wave maximum-

likelihood estimates. This is covered in Chapter IV along with a dis-

cussion of the important properties of the maximum-likelihood estimators.

The fifth chapter briefly outlines some of the considerations to

be made in applying the estimation techniques to electroencephalography.

A limited number of actual dataareanalyzed and the effectiveness of

spectral estimation and filtering is demonstrated.


___















CHAPTER II

MULTIDIMENSIONAL POWER SPECTRAL DENSITY ESTIMATION


Introduction

The multiwave maximum-likelihood estimator presented In this

dissertation requires a planewave filter which effectively reduces

the multiwave estimation problem to a single wave estimation problem.

The planewave filter is realized as a digital filter which is specified

in the frequency-wavenumber domain. The problems encountered for this

type of filter are very similar to those encountered in performing

multidimensional power spectral density estimates. It is thus

convenient to study power spectral estimates before addressing the

filtering problem since the concept of multidimensional space is more

easily explained in terms of power spectra. The primary reason for

developing a method for estimating multidimensional power spectra is

to enable one to quickly search the frequency-wavenumber space for the

presence of planewaves. This then reduces the number of times it is

necessary to perform maximum-likelihood estimates for single waves.

The remainder of this chapter is divided into two parts: one-

dimensional power spectral estimation and multidimensional power

spectral estimation. The first section uses existing techniques for

forming acceptable estimates. Procedures which are extended to the

multidimensional case are particularly emphasized. The multidimensional










estimates are first extended to two dimensions and then to three

dimensions. The actual computer implementation is carried out in three

dimensions. An evaluation of its performance is also made. An
extension to n dimensions is an obvious extension of the work contained

in this chapter.


One-dimensional Spectra
In this section common methods of estimating one-dimensional

power spectral density (PSD) functions are discussed. The important

estimation error considerations are also reviewed. One-dimensional

estimates are used to form the cross-power spectral matrices required
for the maximum-likelihood estimates of Chapter IV. The material gives

the necessary background for the multidimensional case which follows

immediately.

The Fourier transform pair for continuous signals is defined in

the usual way


X(f) = / x(t)e-12ft dt (2-1)


x(t) = I X(f)e+i2nft df (2-2)


for -- < f,t < +- and i=v/T. The units of t are usually understood to
be time (sec) and f to be frequency (sec-1) or Hz. It will always be

assumed that these integrals exist. The functions x(t) and X(f) may be
either deterministic or random. The random processes which are to be

transformed are assumed to be ergodic and, thus, at least wide-sense

stationary.










The power spectral density function can be defined by

+00
P(f) F{R )) S / Rx(r)e-i2ft d (2-3)


where the autocorrelation Rx(T) is given by


Rx(T) = ECx(t) x(t+T)] (2-4)


and E denotes expectation or ensemble average. Power spectral estimates
which first require an estimate of the autocorrelation function and then
an estimate of its transform, as indicated by (2-3), are called indirect
estimates. These estimates are not considered here.
The power spectral density function can be equivalently defined
as [73


P(f) a lim EjXT(f) 2 (2-5)
T-+ T

where |'* denotes absolute value and XT(f) represents a truncated
transform

+T/2 i2
XT(f) A / x(t)ei t dt (2-6)
-T/2

If x(t) and hence X(f) are deterministic, then the operator E may be
omitted and, if the time function is zero for Itl > T, (2-5) becomes


P(f) = 1/T IXT(f)|2


(2-7)










It is important to note that (2-7) is identical to (2-5) if the

power spectral density of a truncated deterministic function is

desired. If, on the other hand, an estimate of the true power spectral

density as given in (2-5) is desired, (2-7) can be used as an estimator.

The resultant estimate will depend on T due to Gibb's phenomenon and,

as will be explained subsequently, techniques are available which

modify (2-7) to reduce the effects of the record length T. Turning now

to the case where x(t) and X(f) are sample functions of random

processes, the same problem exists in truncating the data, but, in

addition, it is now necessary to approximate the ensemble average over

the infinite set {xT(t)} of sample functions.

How these problems are solved to form appropriate estimates of the

power spectral density will now be discussed. Any estimate formed from

the magnitude-squared of the transform, as opposed to the transform of

the autocorrelation function, is referred to as a direct estimate.

If x(t) is a random variable, any estimate of its power spectral

density will also be a random variable. With the estimate of the true

power spectral density, P(f), denoted as P(f), the fractional error or

normalized mean square error is given as [8]


e = E{EP EEP]]2}/E2 [P (2-8)


For the case of White noise (P(f)=constant) E=/Z, which says that the

fractional error is greater than unity. In [9] the complete

calculation of e is carried through for any function with a Gaussian

distribution. The results show that P(O) is the worst estimate with


~








10
respect to f. For values of f not too near zero, the fractional error

E is approximately unity regardless of the length of T. This result

should not be too surprising for, in effect, the average required by

(2-5) is approximated by using one sample function. Estimates formed

in this manner are frequently referred to as raw estimates.

The next logical step is to examine the effects of using more

sample functions to approximate the ensemble average. For comparison

purposes the total record length T is maintained, but the record is

divided into M contiguous segments. This operation is valid because

we have assumed ergodicity and independence of (uncorrelated) segments

provided the segments are sufficiently long. The raw power spectral

density estimate is formed for each segment and the M results are

averaged to form an.estimate of P(f). In general one would expect c

to decrease with increasing M for a given frequency. For the Gaussian

case the approximation derived in [10] is


e2 = 1/M .(2-9)


For the non-Gaussian variates, similar behavior seems probable and

this approximation will be accepted here. If other distributions are

encountered, it should be possible.to derive appropriate error

functions.

The final step in considering power spectral density estimates of

continuous functions is to realize that segmenting a given record of

fixed length, while increasing the stability of the estimate, actually

decreases the frequency resolution of the resultant power spectral








11
density estimate. As is shown in [9] the averaging method described

above is exactly analogous to premultiplying the data record of length

T by a triangular window, Figure 2.1(a), which in turn is exactly

equivalent to convolving (smoothing) the raw estimate by the Fourier

transform of the window function, Figure 2.1(b). For comparison and

later reference, a rectangular window and its transform are shown in

Figure 2.2.

The question naturally arises at this point regarding the

suitability of using alternate functions for increasing statistical

stability. Numerous alternatives have been suggested in the literature

[11,12,13]. One alternative (used in Chapter IV) is called the Hanning

window which, in the sampled data case, results in smoothing with

weights (,4,). The basic objective of any window is to effect a

compromise between improving statistical reliability and maintaining

acceptable resolution. To assist in explaining the concepts of error

control and resolution, the power spectral density of an ideal

sinusoid (extending to -) will be estimated. There is no requirement

for averaging since this function is deterministic, but the effect of

truncation must be considered. The ideal sinusoid and its transform

are shown in Figure 2.3. To limit the time function to T seconds, the

ideal sinusoid is multiplied by the rectangular window of Figure 2.2.

The spectrum of the result, shown in Figure 2.4(b), is the convolution

of the delta functions of Figure 2.3(b) with-the transform of the

rectangle, Figure 2.2(b). The power density is then the square of the

function of Figure 2.4(b) divided by the record length T. From this

figure it should be clear that estimates at very low frequencies tend










to be poorer than at higher frequencies, for fixed T, because of

overlapping sidelobes. This phenomenon is easily controlled by

insuring that the record contains a sufficient number of cycles of the

lowest frequency terms so that the tails of the spectral distributions

are nearly zero at the origin. The second observation is that if a

second time function were added to the first, such that the main lobes

of the spectra of the individual signals overlapped, it would be

practically impossible to resolve two waves from the spectrum of the

sum. This is, of course, the resolution problem. The errors in spectral

analysis arise from the presence of the sidelobes of the main lobes.

They can be misinterpreted as low amplitude signals, or mask the pres-

ence of low energy sinusoids, in the spectrum. This is especially

evident when the effects are shown on a log scale as shown in Figures

2.1(c) and 2.2(c). Ideally, one would alter the raw estimate to

decrease the sidelobes and simultaneously increase the sharpness of the

main lobe. Unfortunately, these demands tend to be conflicting. As a

possible compromise, consider the use of the triangular truncating

window of Figure 2.1 instead of the rectangular window. The new

estimate is simply given as the convolution of Figure 2.1(b) with the

delta functions of Figure 2.3(b). The results are given in Figure

2.5(b). Comparing the estimate shown in Figure 2.5(b) with the former

reveals that the "leakage" of energy from the main lobe into the side-

lobes has been significantly diminished, but the frequency resolution

is only about half the former resolution due to a wider and smaller

main lobe. The optimum choice is a subject of current research [14].


e = .1/TAf (2-10)




































o 0
C

















5--


i 5


c









4-





r-





0 4

C



*g-
... a




t-r

C5




5-

C.
5-4--









C C -
0e
5-




rs CM





5--
U N,



LL_


c &
L a)

CF
04-


0 -24


r ?
C O C*-



*-
C 4.'> 0
4U 0.




4-10
4J
-'4-
's0
Oi






Lo












r- 0
4 r-
O)
mo














Ul "
LO
I 4, -

} Q-
00





0 r-
'-
C

0.. 0
5->




-- + e



CU



Li


-4
















0 0
(a) (b)

Figure 2.3 Fourier transform pair for ideal sinusoid. (a) time
function and (b) frequency function.


(a)


4J





S (b) (fo


Figure 2.4 Fourier transform pair for rectangular-truncated sinusoid.
(a) time function and (b) frequency function.


1.




-fo (
(b)


(a)
(a)


Figure 2.5 Fourier transform pair for triangle-truncated sinusoid.
(a) time function and (b) frequency function.


_ _









where Af is the frequency resolution given as the effective width of

the principal lobe of the frequency domain smoothing function. This

is the equation relating resolution, record length, and error most

generally used and is frequently referred to as the uncertainty

relation [10].

If we now consider the basic time function to be either the sum

of an unknown number of sinusoids or a random variable, it is not

difficult to visualize why the triangular window is generally regarded

as yielding better estimates of the power density than the rectangular

window. In general, any smooth truncating function will tend to

decrease errors due to sidelobes at the expense of diminished

resolution. It is important to realize, however, that raw estimates

due to the rectangular window are acceptable if the true spectrum is

characterized by large, well-separated peaks. In this regard, it is of

interest to note that it is possible to accentuate the main lobe by

actually accentuating the discontinuity (e.g., by differentiating the

rectangular window); however, if this is done, sidelobes can become as

great as the main lobe and the resulting ambiguity seriously limits the

usefulness of such estimates. This is the same problem faced in the

design of weighted delay and sum point array antennas where the

resolution is in bearing angle rather than in frequency [2,15].

The results presented thus far have been justified for continuous
functions only. As this dissertation is concerned exclusively with

sampled data, we shall now discuss the effects of sampling. The

transition is made quite easily because it is well known that sampled

data Ffer a complete representation of a function if the function is










sampled at at least twice the sampling frequency (i.e., at the Nyquist

rate) and acceptable spectral estimates are based upon forming

acceptable approximations of Fourier transforms by using truncated

data. To extend the results to form estimates of power spectral

density functions by using sampled data, the following steps are

explained. The data are assumed to be sampled at at least the Nyquist

rate yielding an N-point time series of length T seconds and time

interval AT. The truncated Fourier transform has now become the

discrete Fourier transform (DFT) defined as


N-i
X(kAF) = 1/N N x(kAT)e-i2ek'/N k=,1, ..., N-1 (2-11)
n=0

where the frequency interval is


Af = 1/T = I/NAT


and conversely


AT = I/NAF = 1/F


where F is the sampling rate. It should now be noted that the discrete

frequency domain function is also sampled. Mathematically this is

equivalent to performing a periodic extension of the original truncated

time series so that the new series is periodic with period N. The DFT

calculates only the first period. Having computed the transform, a

discrete convolution operation may be employed to reduce sidelobes.










The frequency series is thus squared and scaled to obtain the discrete

power spectral estimate.

Having shown that sampled data representation loses none of the
necessary information, it remains only to explain the "rectangular

integration" scheme implied by the DFT. It is recognized that more

sophisticated integration (e.g., trapezoidal or Simpson's Rule) could

be employed to reduce the integration estimation error, but it must be

remembered that these techniques are significantly superior only when

the integrand has relatively sharp peaks and valleys (i.e., contain

high frequency terms). The minimum acceptable sampling, however,

insures that the sampling interval will increase as the frequency

components increase. Increasing the sampling, in turn, increases the

accuracy of the estimated integral. We thus have some assurance that

integration estimation errors are not too sensitive to ill-conditioned

data. We make the usual assumptions that these errors are negligible

in comparison with the truncation and averaging errors mentioned

earlier.

The principal reason for using the DFT in approximating a truncated

Fourier transform stems from the fact that it is possible to utilize

the fast Fourier transform (FFT) algorithm to perform the computation.

This numerical technique was first made popular by [163, is discussed
thoroughly in [17], and a particularly readable account is given in

[183. If (2-7) is computed directly, it is seen that N2 operations are

required (multiply/add). The successive doubling algorithm [17] upon

which the FFT is based permits one to take two Fourier analyses of N/2

points each and combine them in N/2 operations to obtain an analysis








18

of an N-point record of the same data. Successive application of this
technique yields an N-point DFT in approximately N log2 N operations.

A comparison of the execution times to perform a DFT is shown in

Figure 2.6 where it is assumed that execution time is directly
proportional to the number of operations. Graphs for actual computer

execution time can be found in [9]. This great difference:in execution

speed is widely recognized as being responsible for making discrete

Fourier analysis feasible and far outweighs any disadvantage incurred
in not employing other numerical integration techniques. The increased
speed is also responsible for the present trend of forming direct,
rather than indirect, estimates.


T
OFT

107
106

105
104

103
102 F
10

02 1 N
10 102 103


Figure 2.6 Relative execution times.










This section has described various methods of estimating the one-

dimensional power spectral density function. The results of actual

computer implementation are contained in Chapter IV. Details of the
calculations can be found in [19]. The remainder of this chapter is

concerned with extending the concepts to the multidimensional case and,

in so doing, presents a new method of estimating frequency-wavenumber
spectra.


Frequency-wavenumber Power Spectral Density Estimation

The function s(t,r) is a general spatio-temporal function of a

single variable time t and the one-, two- or three-dimensional spatial

vector r = (rx,ryrz). It gives the value of s at time t and vector

position r with respect to the spatio-temporal origin. This function

and its related functions are special cases of general multidimensional

functions. The distinction between multidimensional and spatio-
temporal functions only becomes necessary when s(t,r) is a planewave;

otherwise, s(t,r) may be viewed as a function of an arbitrary number
of variables.

The multidimensional autocorrelation function R(T,p) of s(t,7) is

given as


R(T,p) = E[s(t,r)s(t+r,r+p)] (2-12)


where E denotes the average value, or expectation. Spatial and

temporal stationarity and ergodicity are assumed.

The multidimensional power spectral density function P(f,t) can
be defined as










P(f,t) F(R(T,r)} = / R(Tr)e-12(f+t f) dr d? (2-13)


where the multidimensional autocorrelation function is given by


R(T,r) a E[s(t,r0)s(t+,r0+r)3 (2-14)


For physical systems the spatial vector r is one-, two- or three-
dimensional, T = (r ,ry,rz). The meaning of the wavenumber vector t
will become clear subsequently.
The spatial Fourier transform of P(f,$) yields


i P(f,)e-12w(*) di = Px(f,r) (2-15)
and its inverse is


/ Px(f,)e+i2T(_k ) di = P(f,i) .(2-16)


For spatially discrete two-dimensional systems, r may be interpreted
as the separation between the m-th and n-th sensors, r = (x -x ,y -y ),
and (2-15) becomes

K p+i2T(r )
P(f,1) = z Px(f,x-xny -y)e) (2-17)
m,n=1

where f is the frequency, or,more formally, the temporal frequency,
and I is the "spatial frequency" or the vector wavenumber. f has units
of time-1 and t, length-1. The interpretation of t is made subsequently.









The inverse space Fourier transform of P(f,) yields the cross-
power spectrum, Px(f,p), between two points in space separated by p.


Px(fp) = P(f,t)e+i2(1pj dt (2-18)


The complementary transform may be used to define P(f,1)


P(f,t) = f Px(f,)e-12"1P) d (2-19)


The inverse temporal Fourier transform of Px(f,f) yields the
multidimensional autocorrelation function


R(,) = / Px(f,)e+i2fr df (2-20)


which, when expressed in this manner, demonstrates that R(Tr,) may
also be interpreted as cross-correlation, or covariance, between
points spatially separated by p.
The multidimensional spectrum may be defined directly in a
manner similar to (2-5)


P(f,t) a lim [ EIST(f,k)12/TtTr2n 3 (2-21)
Tt ,Tr

where ST is the truncated n-dimensional Fourier transform of s(t,i).
In order to develop an appreciation of the meaning of the frequency-
wavenumber space, and to introduce the spectral estimation problem,
the two-dimensional case is considered before presenting the theory and
implementation for three dimensions.










Two-dimensional Spectra

Consider a spatio-temporal function of two independent variables,

time and space. Further, let it be required that this function s(t,x)

represents an ideal traveling planewave, which propagates in the

positive x direction. The wave is ideal in the sense that it is not

attenuated or filtered or otherwise contaminated as it travels. The

functional form is the same at each point in space except for the

delay due to a finite propagation velocity. Thus if s(t) is the

functional form at the origin, the value at any arbitrary x is given

by s(t-x/v) where v is the velocity of propagation of the wave in the

positive x direction. For convenience the propagation is usually

expressed in terms of inverse velocity


a = /v|2 (2-22)


Note that is oriented in the direction of propagation of the wave and

represents the vector of delays incurred per unit distance as measured

along each coordinate axis. It is sometimes referred to as the slow-

ness vector. If there is a sinusoidal component of the signal of

frequency f then its vector wavenumber is given by kgfa. The use of

vectors instead of scalers is an obvious extension to the general

three-dimensional space, i.e.,


; (vxVy,vz)
a 4 (ot.ayla
and

t A (kx,ky,kz) (2-23)









If s(t,x) is a monochromatic planewave with inverse velocity a0 and
frequency f0


s(t,x) = A cos [2nf0(t-aOx) + el = A cos [2n(f0t-kox) + 3e (2-24)


where a is uniformly distributed on [0,2]l. The autocorrelation of
s(t,x) is


R(r,px) = E[s(t,x)s(t+T,x+px)3


= A2/2 cos (2nf0T+kOx) (2-25)


From [20] the frequency-wavenumber power spectral density function is
the two-dimensional transform of the autocorrelation

+
P(f,kx) = f/ A2/2 cos (2nf0 +kx0x) dt dx (2-26)


P(f,kx) = 1/2W2 6(f-fo,kx-kxo) (2-27)


The spatio-temporal function is illustrated in Figure 2.7 and its
frequency-wavenumber power spectral density is illustrated in Figure
2.8. Obviously the frequency-wavenumber power spectrum concisely
summarizes the defining parameters of the given planewave. It is
important to realize that if s(t,x) were wideband, instead of mono-
chromatic, the energy distribution in the frequency-wavenumber domain
would be contained in a plane as suggested by Figure 2.9. It is also









24







o





e-
0










to
I i















o
I \i i












\ I












t--
5'-
I I o
C
















I \, Ii









/ 1 '
I I EU
S I 0 o





















^V h x 4.







.- a,

0)





















4-
o









0 Q.

4- i
o


U
C l*





E








S0c



e 8




'4-
= C









EL















94--
/-^ T-
^-^^ *o-
^ a 10
.^' s .
.^^ ^ ^

^^ o A


^~ ^\.01
X '*^ 0- 0
^ '^. u>
-s. <
U- a, t
*- ^s.0
Q. '^^^ j3




^^^ St.
X^
r 10 C
^s. ?
^. 5'?
























4-





IIL.
o 8
X E
> U,

if m

- 0.

0
-D




0.





r 0
a)
"D






I C









S-
t0






u.

a,
In















I,
. 3.







01



U.








27

of interest to note that,if we take the limiting operation, lg4 kO=0
and the familiar one-dimensional PSD is obtained. In practical terms

this condition implies that the velocity is so large that the wave

appears identical at every point in space and that the power spectrum
analysis reduces to an analysis of a function of one variable, time.

If s(t,x) were a sum of wideband planewaves, each with a unique
velocity, each planewave would appear in the spectrum as a power
density function centered on the ray defined by the velocity of the
wave as shown in Figure 2.9. In general the frequency-wavenumber
power spectral density function accounts for the distribution of energy
of an arbitrary spatio-temporal function in terms of frequency and
velocity (or wavenumber).

In considering the spectral estimation problem we again turn to
the direct definition of the frequency-wavenumber power spectral
density function


P(f,kx) = lim EEIST(fkx)I2/4TtTr] (2-28)
TtTr+

where

+Tt/2 +Tr/2
S(f,kx) = i/ s(t,x)e'i2(ft+xx) dt dx (2-29)
-Tt/2 -Tr/2


The problems encountered in forming acceptable estimates for the
two-dimensional case follow directly from the one-dimensional dis-
cussion. For deterministic data the E operator may be ignored;

however, the effects of truncation are present in the form of








28

two-dimensional sidelobes, surrounding the main lobe. If the truncation

window has circular symmetry, or otherwise has simple geometry, the two-

dimensional transform study of Papoulis [213 or Goodman [22] may be

consulted to determine the trade-offs between main lobe resolution and

sidelobe level. These matters are studied in image processing, optical

systems, and two-dimensional antenna theory. However, the lack of

symmetry in spatio-temporal functions makes it difficult to apply these

theoretical results to frequency-wavenumber power estimation. If

s(t,x) is a sample function of a random process, one will usually want

to know the relationship between the number of spatio-temporal samples

needed to approximate the ensemble average, the window one uses for

truncation and, of course, the statistical error function. While these

matters are quite important in determining the value of a general power

density estimator, they are not investigated here for two reasons.

First the power spectral density estimator required for this study is

needed only to assist in the detection and estimation of strong signals.

If spatio-temporal fields are encountered which have fairly flat

frequency-wavenumber distribution, the power spectral density estimator

simply is not used. The second reason for not surveying the statistical

advantage of averaging and/or windowing stems from the fact that

frequent applications are characterized by arrays of very limited

spatial population. The array considered in Chapter V consists of a

4 x 4 array. Thus if an attempt is made to smooth using existing

techniques, the reduced resolution renders the estimate virtually

useless.










Three-dimensional Spectra

The spectrum analysis of spatio-temporal functions of one spatial

variable and time is conceptually easily extended to the general case

of three spatial variables plus time. In the latter case the wave-

number k becomes the vector wavenumber I which has three coordinate

components, i.e.,


PSD{X(t,r)} = P(f,t) (2-30)

where

S= (rx,ryrz}
and

S= {kx,ky,k (2-31)


From a practical viewpoint one spatial coordinate is too restrictive

to be of much use in modeling physical systems and at the other extreme,

with three spatial coordinates, an element in the frequency-wavenumber

space is an element in a five-dimensional space whose coordinates are

power density, frequency, kx, ky and kz. The data processing task for

handling data of this magnitude is considerable as is the task of

interpreting the results. As a compromise it is common in sonar and
seismology to use two spatial coordinates. This compromise is also used

in Chapter V. An examination of this estimation problem follows.

Consider a function similar to (2-24) but extended to two spatial

coordinates


s(t,x,y) = A cos (2nf0(t-ax0x + aygy) + e)


(2-32)










or using vector notation


s(t,-) = A cos (2nf0(t-ao.r) +e)


s(t,Jr = A cos (2nft 10.r + e)


where t0 = fa. The
delta function


K 6(f-fo0, `)

The vector velocity


frequency-wavenumber spectrum of this becomes a






and inverse phase velocity are related by


The term phase is added to emphasize the fact that we are dealing with
a projection of the actual (group) velocity onto the x-y plane.
In attempting to sketch the frequency-wavenumber spectrum, we will
use points to indicate three-dimensional delta functions, thus the
example cited may appear as in Figure 2.10. Note that the direction of
propagation of the wave is the same as t0 if the kx axis is relabled x
and the ky axis is relabled y. The magnitude velocity of the wave is
given by the slope of the line from the origin to the dot, i.e.,


(2-34)


Iv0l = 1/|a(0 = f0/kO


(2-33)


* -).I -"I
a = V/|v|


1 =f6,




























o
-S-



4- 0-' a

o 4 0

41 w n
0u t CLr-




0 0




r- 4in




E"
OL. 0





CD
E u


r 0














CL
'~ C-







J3
4-


4-





1 I
ci



\ '*0 .
'f.--- -.-------------------------- ---------Co0





\ r 3-







LL. 4





O C.
4-U



\..
\ / '-










Obviously then the power spectrum characterizes the vector velocity,

frequency, and amplitude of the monochromatic wave. If wideband waves

are considered,the points must spread along the line passing from the

origin to the dot.

To enhance our knowledge of this case we compare our technique to

that used by others. As noted in [233 beam-forming and one-dimensional

spectral analysis can be combined to accomplish a study of the velocity

and frequency structure of the wave without any direct suggestion that

estimates of frequency-wavenumber spectra are being obtained. This is

referred to as the time domain method and is illustrated in Figure 2.11.

The spatio-temporal function is sampled as indicated by the inputs si.

Each channel is delayed to align any incoming wave with a specified

vector velocity. The channels are then summed to diminish the effects

of noise. The resulting function is referred to as a beam and is then

spectral analyzed in a convenient manner and the end result is power

as a function of (foo0). While the method is conceptually simple, it

is subject to severe limitations in velocity filtering or weighted

delay and sum operations. Generally, the weighting is adjusted to

suppress interfering events when the frequency composition and bearing

of the unwanted events are known a priori. With no a priori statistics

maximum-likelihood weights can be obtained but, in practice, this offers

little advantage over straight delay and sum with unit weights 13].

The next method, referred to [23] as the beam-forming method,

frequency domain, is termed the conventional method by Capon [3]. It

consists of an extension of (2-18), i.e.,
















































In



I 0m
I W




L--


~~~~~.1


tL +S- +s-
-4J +4 ... 4

(n l vv)












P(fst) = Wp Px(f,)e-i2DP (2-35)
p

where Px is an estimate of the cross-power spectral density and W is
a window used to improve the estimate of the two-dimensional transform

of Px. The advantage of this technique over the former is that it is
not necessary to form a large number of beams. From an implementation

viewpoint this estimate is appealing as standard techniques are available
for estimating Px and the two-dimensional transform can be estimated

using windows which may be symmetrical as the units in both dimensions

are identical. The estimation errors associated with this procedure

are covered in [3,233.

Whenever estimates involve approximating Fourier transforms, one

encounters the problem of imperfect (not ideal) resolution in the

transform domain. The usual techniques involve choosing a truncation

window in one domain which yields acceptable resolution in the transform

domain. Capon [20] developed the "high resolution" estimate which

utilizes decision theoretic ideas to form an optimum estimate. By
definition the high resolution is


P'(f,k) E [ q q (f)ei2 k(j-x ] -1 (2-36)
j,t=1

where q is the inverse of the estimated cross-power matrix. The essence
of this technique will be covered in Chapter IV along with maximum-

likelihood estimation. It provides excellent estimates, in general, but

is computationally inefficient.










The estimation presented in this study is based upon estimating
directly from a form of the defining equation


P(f,) lim 1/16TtTxTy El| s(t,r)ei2I(ft+k'r) dt dil2 (2-37)


where t = TtxT y. As an estimate this becomes


P(f,t) = 1/161 IDFT(s(t,')})2 (2-38)


This method of estimating the frequency-wavenumber spectrum will be
referred to as the direct FFT method. This procedure is similar to that
advanced by Smart and Flinn [24]. The difference is mainly due to the
method of calculation. The reason for choosing this estimate is to
utilize the computational efficiency of the multidimensional fast
Fourier transform [16,173. While reference to the multidimensional FFT
is contained in the landmark paper of FFT [9], the literature seems
devoid of its application to frequency-wavenumber spectrum estimation.
This may possibly be explained by the requirement of the FFT to utilize
uniform sampling in each coordinate direction. It should be noted,
however, that sampling intervals in each direction need not be equal.
It is also interesting to note that early attempts at one-
dimensional power estimates were based on the indirect method of first
estimating an intermediate function covariancee) and then transforming
to obtain an estimate of power density [9]. With the introduction of
the FFT, the indirect method was replaced by the direct method which
essentially estimat'- the magnitude-squared of the Fourier transform.










Where the requirement for uniform sampling can be accepted, the use of

the direct method of estimating frequency-wavenumber power spectral

density via the multidimensional FFT will very probably become

increasingly popular.

As mentioned previously, the question of finding suitable three-

dimensional windows to improve the raw estimates is not investigated

here. The Fortran programs written to implement the estimate (2-38),

however, include provisions for multidimensional windows. This point

will be covered in the next section. Whether or not special windows

are included, the computations required by direct frequency-wavenumber

PSD estimation are significantly less than those required by other

estimates. The frequency domain method appears to offer the most

efficient second choice. In order to compare these two estimates,

assume we have N1 x N2 x N3 elements of a spatio-temporal function to

be analyzed, where the ordering is time by x by y. The direct FFT

method can be accomplished in about NIN2N3log2NIN2N3operations [17].

Assuming the one-dimensional FFT is used, the frequency domain method

requires N1log2N1 operations to transform the N1 time functions plus

(N2N3)2N1 operations to form the cross-power matrix and finally N 2 N

operations to perform the final two-dimensional transform. Table 1

compares the number of operations required for the two methods for four

particular values of N1, N2, and N3. The execution time for the

32 x 4 x 4 case on an IBM 360/65 computer is about 20 seconds using
the direct FFT method.











Table 1

Number of operations

N1,N2,N3 Frequency Direct
domain FFT
method method


32,4,4 8.6 K 4.6 K

32,8,8 135 K 22 K

1024,4,4 266 K 65 K

1024,8,8 4100 K 130 K






Implementation

Sets of Fortran subroutines have been prepared for the purpose of

estimating the three-dimensional frequency-wavenumber power spectral

density functions via the direct FFT method. The dimension and format

specifications are specifically written to handle input arrays of size

32 x 4 x 4 in the order of time by x direction by y direction. This

corresponds to a typical array used in other studies to which this

research has been applied. Straightforward programming modification

could be incorporated to handle virtually any array size.

The explanation of the purpose of the principal programming units

used to analyze simulated data follows. The analysis of actual data is

given in Chapter V. Actual programs are included as Appendix D. The

flow chart of Figure 2.12 should be noted in the following discussion.











Define wave as time function, select
Call PWAVE I...... vector velocity, and select signal-to-
I noise ratio (SNR).


Call PLOT16 ...... No options. Utility routine, plots
all sixteen channels on a single page.


Call NEGATE ...... No options. Rotates function about the
x and y axes.


Call FWNO ...... Select time domain window or
frequency domain smoothing.


Select frequency range and interval
Call PLT3 ...... for plotted output, select form of
output: printed number array, grey
level array or both.


Call NEGATE ...... This call nullifies the effect of
the first Call NEGATE.


Figure 2.12 Simplified flow chartfor simulated data.


Subroutine PWAVE accepts a given function of time which is taken

to be the signal at the origin in space. Fifteen additional functions

are then generated which correspond to the time functions which would

be received by an array of sensors for the specified vector velocity.

The array is assumed to be square with 2 cm sensor spacing. If.ambient

channel noise is desired the signal-to-noise ratio (SNR) is set to

0 s SNR 990. If this condition is set each channel receives additive










independent Gaussian noise [19] with a peak rms noise of SNR. If

SNR > 990 no noise is added. If SNR s 0 noise only is returned with

unity rms value. PWAVE may be called repeatedly to generate a sum

of several planewaves.

If we refer again to Figure 2.12, PLOT16 is next to be called.

This routine provides a quick look at all sixteen channels with

relatively low amplitude resolution.

Subroutine NEGATE is used to transfer s(t,x,y) to s(t,-x,-y).

This is required whenever a standard FFT is used in which the exponential

term of the forward transform is given as


e-i2n(ft+kxx+kyy)



In order for the t and v to have the same direction, the exponential
term must be


e-i27(ft-kxx-kyy)



which is effectively accomplished by NEGATE. If NEGATE is omitted

the kx and k directions will be reversed in the final output.

Subroutine FWNO estimates the three-dimensional spectrum basically

by calculating the magnitude-squared of the FFT of the input spatio-

temporal series. To enhance statistical reliability and reduce

leakage several options are available. One may specify tapering

(windowing) in the spatio-temporal domain or smoothing of the raw

spectrum by either Hanning or Hamming [9] smoothing. The frequency








40
domain smoothing only smooths in the temporal frequency direction

and not in the wavenumber coordinates. The tapering utilizes a

10 percent cosine window in the time coordinates and the spatial

planes are bordered with 0.5 weighting with the center four

elements weighted at 1.0. FWNO may be called repeatedly, with

different samples of the same process and the results averaged.

To increase resolution in the wavenumber domain zeros are

added in the spatial domain. Each x-y spatial plane is padded

with zeros to double the resolution in the transform domain.

As adequate resolution is available in the temporal frequency

coordinate, zeros are not added to the time coordinate. The

next step is to calculate raw power estimates for positive

temporal frequencies. IBM's subroutine HARM is used for the

actual FFT [25]. As the ordering of the spatial and temporal

frequencies from any FFT is positive frequencies followed by

negative frequencies in all coordinates, the last step reorders

the frequencies to normal ascending order.

If the input series is a spatio-temporal impulse response,

the above applies except that the zeros are added to the spa-

tial tails of the impulse response. Instead of bordering the

x-y plane with zeros, zeros are inserted at the positive and

negative ends of the spatial records. This optional method is

selected by calling FWN05, instead of FWNO. FWNO is actually

an entry into FWN05.

Subroutine FWNO returns a power density array which may be

used for any purpose. If a line printer display of the data is

desired, subroutine PLT3 can be used. This routine outputs








41
power as a function of kx and k for each temporal frequency

desired. The output is normalized on a negative decibel scale

with respect to the maximum of the input power array; At each

temporal frequency the output consists of a matrix of numbers

(db's) whose position corresponds to the proper position in

the kx-ky space. Contour plots, using grey levels as described

in [26,273, can also be added. The contour plot more than

doubles the cost of the output, but renders the results much

easier to interpret. If normalization is desired with respect

to a previous PLT3 call, PLT4 should be called. This is useful

for observing the effects of filtering, comparing noise level

to signal level, etc. If normalization is desired with respect

to a particular point of the input array, PLT5 should be called.

This is useful in examining the amplitude-squared response of

a three-dimensional filter as it allows one to normalize with

respect to a specific point in the passband. Ripples are then

seen as positive or negative deviations from 0 db.

If it is required to return the spatio-temporal series to

its original state, NEGATE must be called again. NEGATE was

not included within FWNO as several calls may be made to FWNO

without ever returning to the spatio-temporal domain.


Results

The usefulness of frequency-wavenumber power spectral

estimates is demonstrated by application throughout this

study. In this section, an example of the application to the

analysis of simulated data is given. A general assessment of








42
the performance of the estimate is also included.

To illustrate the performance of the Fortran subroutines,

seven truncated, monochromatic planewaves were generated and

summed using PWAVE. The frequency-wavenumber spectrum of the

resulting spatio-temporal function was estimated by FWNO and

plotted via PLT3. The input functions are listed in Table 2.

The temporal sampling rate was 50 Hz, which yields a

folding frequency of 25 Hz and a time interval of 0.02 sec.

The record length was 32 points or 640 msec. Spatial sampling

took the form of a 4 x 4 square array with a spatial interval

of 2 cm. The "spatial folding frequency" was then 0.25 cm'1

and wavenumber interval, 0.125 cm-1. The spatial dimensions

of the array are 6 cm x 6 cm.

The raw spectral estimates were computed by FWNO and the

result listed by PLT3. The actual output is reproduced in

Figures 2.13 through 2.29, inclusive. Each figure corresponds

to a segment of the frequency-wavenumber space for which the

temporal frequency is constant. Thus Figure 2.13 gives power

as a function of kx and k for 0 Hz. The power level is re-

corded at each point in the space in negative db with respect

to recorded peak power density. The minus signs are suppressed

for simplicity. The peak power density is calculated for an

input given in micro-volts or micro-amperes. The power density

scale is limited at -90 db. Power levels below this point are

recorded as -90 db. The coding of the power density regions

into grly tones is given in Table 3. The coding scheme was

chosen to emphasize the higher power level and, in particular,










Table 2

Listing of overlapping planewaves

Number f0 V0 )VO
(Hz) (cm/sec) (degrees)

1 3.1 30.0 120.0

2 7.8 42.0 90.0

3 7.8 125.0 90.0

4 12.5 75.0 45.0

5 12.5 75.0 -45.0

6 17.2 100.0 0.0

7 21.9 100.0 0.0






Table 3

Grey levels

Overstruck
Level db range symbol
symbols


+w, -3

-4, -9
-10,-15

-16,-21

-22,-27

-28,-33

-34,-39
-40,- =


( , )

( , ,-)

( , =)

( ,+,+)
( ,X,X)

( ,x,x,=)

( ,*,X,0)

(#,M,W,O)
























PI HATIRl F(IQ TFMPDOAL FREDO 0.0 HZ
IN -0I8 .9.T. PLAK PODrn AT 7.8 HZ
OF 0.ib1399i 04 PICOWATroSEC*CM*2

30 39 41 36 41 39 3R 45


34 30 37 35 31 36 34 37


3. 29 36 35 29 29 30 37


34 34 36 37 35 29 31 51


37 38 43 30 43 38 3? 45


31 29 35 37 36 34 34 51


30 29 29 35 36 29 32 37


34 36 31 35 37 30 34 37


VELOC


0.250 0.0


0.108 0.0


0.125 0.0


0.063 0.0


0.0 ******1


-0.063 0.0


-0.125 0.0


-0O.I4 0.0


-3 -2 -1 0 1 2 3 4

*ESSSeSSOMIIIIIsIssIseEsaEmastII IIIE



==06amalgamgogeame sMRW a smemo miuea

MES4X5M5N*ea55sag5N5 4MM8SxMaSE a555




mas5 am nun EamanaEESge K::: : 8:::5Boo
0S1SLOCSASLMAN:SSESL0ASSnama0SSl Ell3

MMMMMMMENASAWHASIAAI8BMWAAAAeMIMagII
4585588 ESSESESSEEWSEN 5555555 ASSAA





WNEIIS5I5N 55ExxISBa IXIEM W 185850
ESESNNEMMaSISnNNSSSSUDIEEE'SEANNaSER
*8agsss48SslSESSSEEUSES SUSSIonAMBAm






At LOCAL 4N8: VELOC 0.0 CM/SEC. UESRING n-IlS. OEGREF5


Figure 2.13 FWNO output, f = 0.0 Hz.

























PU MAT1lX FOR TLMPORAL FRO0 1.b HZ
IN -08 m.R.T. PFAK DOW*E AT 7.8 HZ K VELDC
OF 0.1961399E 04 PICONATT*GEC*CNM*2

39 32 32 34 35 30 36 42 4 0*250 6*25


30 27 38 29 2B 37 30 36 3 0.188 8.33


29 30 29 20 29 26 20 39 2 0.125 12.50


33 30 25 34 27 24 3 2 1 0*063 25.00


3? 3 3232 3 27 29 36 34 0 0.0 *4*****


30 32 30 29 *1 32 33 42 -I -0.063 -25.00


31 28 27 36 32 20 35 33 -2 -0.125 -12.50


36 35 3R 42 34 33 34 37 -3 -O*.I8 -0.3J



-3 -2 -1 0 I 2 3 4


c*RKxxx lxxnKmiX iMa lia
**KUXCMiiAiAlAXAMOMAx0UExsmeMSEiSl0l

AT LAMCA MAX MANVCC 11*A CM/ BEARING 27. REES

XAMANNlAMMMAMMAMAMAM AMAu MAIAMMNAS

A AMAKAAM NMAMMA MAMxxMUM KMA NMAMMA
AMMAXAMMXMMAMMASSMMMAlMMNAMMMMMMlNAI

AA MARrASM MAMAMAM AMA KM M*AMMAMMA A





MKM MM AMA MA AM MABMAS8ANAAMKM MANH
AIIIAMMIIMRIAMMEMAMMIXMMMMIAA OXOBMO
RMMERRER1MMX1KKXKXXRXxMBMWAB
MAMAMAMM MM MNMJMMSCIIMIIMMAM MA A800







MMEUSSMMMSISMUMSSSSMSSaAAUMMeMIadI


Figure 2.14 FWNO output, f = 1.6 Hz.













46












PM MATRIX FOR TEMPORAL FPWO 3.1 HZ
IN -08 w*R.T. PEAK DOWER AT 7.8 H K VELOC
OF O*lQ6l39a E 04 PICOMAT TSEC*CN**2

32 21 I* IS 23 20 27 27 4 D0.250 1250


32 5I 12 IT 27 22 25 19 3 0.100 I6.67


23 7 2 4 18 16 34 13 2 0.125 25.00


22 A /2 3 1 16 26 14 1 0*063 50.00


30 16 9 17 24 23 25 0 0.0 *******


34 23 1 22 / 26 27 25 -1 -0,063 -50.00


31 26 18 10 24 33 33 33 -2 -0.126 -*5.00


34 20 16 20 36 25 34 25 -3 -0.I50 -16A67



-3 -2 -1 0 1 2 3 4

HMRKXXMeAmn.ss *XMMNXMMAMXXXXMXMNMM
MNEA **** ::****XXMM*XMMMAMXXXRXRM
MMmm*A*a-..m. **MMMMMMAAMI= MMXXXXIMM
MMNXXr4asaas-......IXXXXNKXXXXXXXXX*

X:t --w-------m----n=-+*******KMM,
.. .M ......AA MOX

AMXt ----- --- *MI XXXXKNXMMES
XNM.M .--- X*M XXXMMMMXXxNMM**







XAXXX *AM*9AAAMxAM:: KXXXAXXXMNA*
MxrKAXK---- --MXXXXKKKB*MMNXX
AT LOCAL MA: VLOC 3----536 CEC EARIN EGRES
MAMAtwo --- ---wAn* AAMANMMMMMM
MA.....99---- t MA DMAMMMMMAMA.
AMME ...senAMM*XXXE MMMMAMMAXA
MMMAAxA* AAA*M*MMMMAMM RxxxxMMMMMMME
EMMAMAMAK *M$MARMMMUS MMMxMMMMMRXXM

MMMMMMMAA**A**A AMMMMMMMEMEMMM
MMROAMM*AXM3569MMERMMERM MMMMG 3MMMM
AM MM M*AM *M AMAMAMMMM MAMMMENMAMM A
A MNNM MMA MM MAMMAM M** A MARARM EMA


AT LOCAL 4MM: VELOC 35.36 CM/SEC. HEARING 135. DEGREES


Figure 2.15 FWNO output, f = 3.1 Hz.


I

























PR MATRIX FUR TEMPORAL FPREI 4.7 HZ
IN -08 4R.*T. PFAK pnwCq AT .8 HZ"z
OF 0.1961399E 04 PICODATT*SEC5CMC*2

20 25 35 25 25 39 30 38


20 20 29 22 24 31 29 36


33 25 24 32 22 22 36 29


33 23 28 22 17 21 31 28


28 20 29 19 20 29 29 38


30 28 26 27 37 31 33 36


34 29 34 32 25 20 30 34


36 33 33 34 38 36 35 45


K VELOC


0.250 10.75


0.18a 25.00


0.125 37.10


0.063 75.00


0*0 000****


-0.063 -75.00


-0.125 -37.50


-0.108 -25.DO


-3 -2 -a 0 1 2 3 4

HMRMS 88 M88 M 8 888888 888 oB MM38833


MMERXMERNMNNNMINM NMN.NNM: ::XMARNRMA

8 8 8888 8 88 8 88 N88 8 x8 33388 8 83

x8xx 8 xx4 xx8 xx8 88 xxxx8*** 8m88Rm3


NXMNU AINKXNX ......*I*rXxx qxxx** mmo
8Xn**IIXXXMXXX#MX6XX6XX****MMMMMHN


383XNNxxmxxxxxxe8sA.**8*xxu8MnM8888
AT8A8 MAN:8X5XX30N.A88888888NMM8888MM
888888888888885888888 8X8888R88X83O36

IKMMNNNXRRXAIONEN*IRMNMEMMKMNMEMMEXN
8NXNRAMMEIaXxxIIxMMKatxaxMMERMO8 2AA
IRNRRNRXaXXXaMNMRNMNMIauxIanX*XMgARM

8i88lxxx88N88N8N888M888MM88M8M888883



AT LOL8 L 8 A8 : VELOC 8 53803 C8/8S C BEn8NG 8 858 OEGE~ES


Figure 2.16 FWNO output, f = 4.7 Hz.












48










PW MATRIX FOR TEMPURAL FRE 8.63 NZ
IN -OB M.R.T. PCAK POWER AT .8 Y3" K VELOC
OF 0.1963199E 04 PICOWATTXSEC*CM**2

26 28 31 22 24 63 30 37 4 0.250 25.00


29 29 25 22 24 30 36 33 3 0.18 33.33


35 22 25 29 20 22 35 30 2 0.12S 50.00


30 24 29 20 17 21 28 31 A 0.063 100.00


29 32 28 20 21 27 30 40 0 0.0 *******


33 30 30 34 37 35 36 *4 -1 -0.063 -100.00


39 36 34 32 26 28 37 39 -2 -0*125 -50*00


20 2? 36 35 38 35 32 39 -3 -0.1s0 -33.33



-3 -2 -1 0 I 2 3 4

4XHXX400XMXKXR8KKMN8XMN3m03330MMMN0

EXX4XKARXXRXKX XXIXXXXOO3UOOUO3030
RMMXxXXXXXNXXXXIXXXXM ..XX X. NM XX M33
KxxXXMXMM xxxxxxxxxsxxzXXXXXMMUMAO
AMXRMXRKXXMRKXEAXXX XXM XM MNMM 4X

X0MXXXXRXXXXRRXXX... XXXX4KMMMM
MMANENxKNMMMMMMANNXMMAM xxMxM nM
MMNMMXXMXxxxxxxx**e****XXMe xxxxNNMW
XXXXXXXXXXNNaXRX4XX. XXXXXXXXX048000
X0MMiXNXxKxXXXXMMXI*****.**XXXX MXXO

XXIIXKNMKlNxxx***** XXXXXXXMMXMxX0m8





AT LOCAL MAX: VELOC 70.71 C8EC* BEARING .5. OREES
0KMMM MXXMM8M XXX XNMX XX MER O ff3 03
gIXRU 3IXR XaX11XXEEnXXMMEXXOOII00IB
44 48 XX08 000 00300 008 MEME 300083
BIBIXMMMIIMKI#NRNMMMENNXXXXNKXNMMMII


E84NXXNR3IMOIIIIIO3IOIIXIIII0MEEEOO


Al LOCAL MAX! VELOC ?O.?I Cs/SEC. HEARING 45. OLGREES


Figure 2.17 FWNO output, f = 6.3 Hz.


























PM MATRIX FOR TEMPORAL FREQ 7.6 nZ
IN -DB w.R.T. PEAK PUWFP AT 7.8 HM K VELOC
OF 01961399E 04 PICODATT4SEC*CM**2

13 32 5 2 6 SO 14 34 4 0.250 31.25


II 25 3 4-- 32 II 36 3 0.180 41.67


s 24 2 14 35 2 0.126 62.50


II 29 20 12 34 I 0.063 125.00


12 34 3 25 12 35 0 0.0 C******


2 32 20 1 20 44 6 4 -0063 -125.00


20 20 II 7 10 2 20 30 -2 -0.125 -62.50


31 28 24 22 J 31 35 33 -3 -0. -4.67



-3 -2 0 I 2 3 4

m N .OX.--. ----...O. N .XOMS

m**eXXXA*- ----X*XMNXXS**XX
*OAAXNftlr -AXNNXSfl4XNU

.***5 0 m1-- --As* X4A.xxNU

S .....---
ANONxa --* S*t4 *440X

S*x.-- 44


**OANMN* -4 4XOMSa *XSS
XXX XXX M#

XX ONMMXXAA .A* SMXUXNmm OXOXEE
X05 NMIZXS ... .=. sAAommN OXXNNIIX

---- 14**------=40 XOXMMMNX
AXIOM XXXOS



AT LOCAL MAX: VELOC 41.67 CM/SEC. HEARING 90. OLDECGREES


Figure 2.18 FWNO output, f = 7.8 Hz.












50










PW MATRIX FOR TEMPORAL FREGQ 9.4 HZ
IN -00 W.R.T. PEAK PDOMER AT 7. HZ K VELOC
OF 0.1961399E 06 PICOWATTOSEC*CM**2

30 28 27 22 24 39 31 39 4 0.250 37.50


27 25 37 28 26 34 34 34 0.166 50.00


30 20 27 21 20 23 27 39 2 0.125 75.00


30 35 27 21 19 20 24 32 1 0.063 50.00


32 33 38 29 23 23 29 32 0 0.0 ******4


33 20 30 31 33 56 37 46 -1 -0.063 -150.00


29 27 36 35 35 28 31 35 2 -0*.25 -75.00


27 33 29 23 25 31 39 33 -3 -0.10 -50.00



-3 -2 -1 0 1 2 3 4


...NX...MMNMMMRNXNXXXxNMX XMNBAMMNBm



MEMNNMMNXI NaMMMMMXXMaNIeXaggaag00

















AT LOCAL MAX: V6LOC 6 106.07 CM6S6C. REARING a45. DEGREES
6666R66R6666666666 6666xx6xx6x66M663
666666666666666666666666666666666666
MMXXMNMMlNMElxx++++++++XXXMXXNNNMER

MNRnI|aaxaxXn*e**XIX**+XXXXXONXXxXxa




XNIana*xxxxxuxxxaxaxx****oo+ooExxonu


NMxx*MMMABxAxxMoMrMMMxxxxxxxxagao


MMMMMMXIMMXNMXIXIXXXKXXI**I*I*NNNtHW

AT LOCAL 6 AX: VEL6C 6 10606 C66SEC6 RE6RING 6 5. O6GRE6S


Figure 2.19 FWNO output, f = 9.4 Hz.
























pN NATRIX FOR TEMPORAL FREO *10.9 HN
IN -DO W.R.T. PEAK POWER AT ?.8 "Z
OF 0.1961399E 04 PIC0WATT*SEC*CM**2

2? 27 33 23 23 32 33 36


20 28 33 25 26 34 30 41


36 36 20 21 21 23 26 34


29 36 28 23 21 20 23 30


35 30 30 40 24 23 26 35


30 27 30 29 29 38 42 39


26 31 44 32 33 30 33 34


S6 36 30 23 23 30 s0 43


K VELDC


0.280 43.75


0.100 56.33


0.125 87.00


0.063 175.00


0.0 *****. *


-0.063 -iT7.C3


-0.125 -07.50


-o.1ea -s.*33


-3 -2 -1 0 1 2 3 4

MNMXKKXIIIXIIKM XXXKMMIMNNXNM KKXrBaBB
IIXWKKIMIIXLXHMXXaXKaKK1KHIIKMMBBBa
MxxxEM l:xlc xxxxxxa lxxuxaa aptBB gAAoeg

xMMJMMxxX MKNKEXKXXKXXX Xl BR HABi OBA


IXBMXMMXM***I****XXKXXXXXNXKMNNMM


MMIIMILMXXXMXXKMMKKX**MMX**XXXXXXXKM


AT LOCAMMA N.VE-A RMMMN XR C EAIN 2. EGEE

asstuseasemeBAnN RntKaaNENMEMNMESSI


MMARttmRmmma1NMMxxx1*x11xxMMEMMMADag
BRA1B1EENRmmg1iMMM xxx N MMRXXXBet




AT LOCAL MA": VYELC = 7.2& C/15EC UlARING = 27. DEGREES


Figure 2.20 FWNO output, f = 10.9 Hz.


























PW MATRIX FOR TEMPORAL FREO .12.5 HZ
IN -OD .P*Rer PFAX POWER AT 7.0 HZ
OF O.IS4I394E 04 PICOEATT*SEC*CME*2

27 30 30 26 25 31 41 31


14 32 14 25 a 4 7 43


12 43 12 22 5 I 5 31


16 33 7 26 5 6 30


35 30 20 39 2S 23 27 49


14 29 I5 32 3 7 34


12 35 s 34 5 4 30


15 37 I 24 a 4 7 34


K VELOC


0.250 50.00


0.3I0 66.e6


0.125 100.00


0.063 200.00


0.0 *******


-0.063 -200.00


-0.125 -100.00


-0.16s -66.67


-3 -2 -I 0 a& 2 3


ENNrXKN EWExx EM NxxJEE NEa- --X-XNE
*..ANKNE AXXEf ... =..-*A..XE.n.

*EUU E*X***tXX"--- -XX.-f*** a*a

songagan --- -----~ag
*eMAW *WEEAEA w-N- --ttx,
aE NEE***B----E -----*Xm
SEEUNEEMN*XNNX**NXE------N---E-NEXN
SENEIEEE*EEXXE=aE-E----- ---E--*X


*ENXNNMMXEEMMMmEMEXEAEa----mNiNE
NN ManxxXXM:MIMImXImX****'*** MX

....MNI..X*XXXI IX++.. ++a* *:X.::
***l lxxt. "XXtXNM*=----- --a XN8
**X MM**:::m ***XMNxX --- ----
*tE ENXEX**EXENNM--- ---ANNBE

*EXMEEXANASEA#XEK+*N--- --N*XNE
TIXXALx=WAX s=----- ---C--SE A
w+XNMD~jta*.+++xa++.--:: --------- X..

AT LOCAL AX: VELEIC 70.71 CMSEC. BEARING -45. DEGREES


Figure 2.21 FWNO output, f = 12.5 Hz.












53











Pt MATRIX FOR TEMPORAL FREO 314.1 HZ
IN -OB wReT. PEAK POMEN AT 7.e HZ VELOC
OF 0.1913990E Oa PICOaATT*SECMCM**2

31 33 33 31 28 30 43 33 4 0.250 56.25


43 36 30 27 27 29 33 39 3 0.183 75.00


36 45 30 25 23 2* 26 31 2 0.125 112.50


34 33 30 30 2b 22 2* 32 I 0.063 225.00


37 31 28 36 26 23 27 B4 0 0.0 ttl****


4? Jl 29 3* 29 28 3 3 30 -1 -0.063 -225.00


36 36 31 3* 30 27 32 35 -2 -0.125 -11t.SO


29 *1 3* 28 2* 25 35 31 -3 -0.18 -7S.00



-3 -2 -1 0 I 2 3 4

aNMtR MMaa RN ara auxx ax ngamasama

*iriinMririiNxxxxxxxNxxxxRx.n.Mamna
ImLIrMKXRXNKNNEXXXaNXKRIXKXXNMMMND


BMxNIaSxNsR**XraM*NRaxMNNaaaxaxalN
555 955 N N KM Eta ta a MNMEN MMEMN


SMMMMMMMMNXXXXXXMXXMXXXXXXXEMXXXXNII

ARM LnCAXL AxxMA VELOC a M 10 2 CSEC. IG 7
55535 MM MRNMMM M NMaxxxxxtxXEMMN MM
MMaNN MtNMNtMMSttIXtXtXX XNN NItMME
EMElA M MN Max aANAmM Aaat aXtNMNMNMM


SMMaMtNRXMMMUMXMEXXRXXXNNXXaMN*MNIUM
IB*HMNHMtUMMXIatMRRNMtEMNMMMEMSNBSII
IEMMMiNaIMXSXMNRIXIXNMMMMMMIINaXN*MM
Rmigat tat atxM eNMM MNMM NaxamM s
IMammammaxxaxaxaaaSAMaNMNMaxaxNMMUU


a awls M MM asxatM aNItMX MMM MM'MMMa
M IMMAISMMMMMa ata X aalXIMaa MOOMNAMM
AMtLICELM A5NRNaE 52RE EAIGX XX27NMDNE E

At LOLAL MAR: VELOC 100.62 CM/5EC. BEARING = 271 DEGREES


Figure 2.22 FWNO output, f = 14.1 Hz.
























PW MATRIK FOR TEMPORAL FREO *IA56 HZ
IN -DB .<9.T. PEAM PONWF AT 7.0 HZ
OF 0.196l399E 04 n COIATTSEC*COC*2

3S 36 32 37 35 30 32 34


0I J3 32 30 30 29 29 34


3? 52 34 29 26 24 26 33


3J 36 30 34 28 23 25 34


&1 32 27 33 29 24 27 39


36 33 26 30 33 2? 33 39


31 38 29 32 34 28 33 35


30 45 33 35 29 26 31 34


K VELOC


0.250 62.50


0.380 83.33


0.125 15*.00


0.063 250.00


0.0 00*****


-0.063 -250.00


-0.125 -125.00


-0.188 -83.33


-3 -2 -1 0 1 2 3 4

0e6sse*066M56g603B60A6R 60N0MN6NNN00
AGSfJIBIDXMNWWmWMN ANNNMMANNXN NNKNN

BIWIAAIMIIMEMMMMMAOIMMMANRXNXI*RNNI
iBAIMMNMMMAMMMMMMXMMNXXA AMXXMA




WatWMANMMMWNAN*X*MMII1axivxMMA.:a
IIIIIIIBEIOIXXMEMMMMMIIIIXMXMERRNIRH





2a=2exNxxxxxxxxxxxngxaxxxxxxx: ana
BAMMMMMMEMRREXNXMNMKNMAMMMMKNNOBAS



MOMAIIItXX*II*MEMMXENXI*NXXNNNNERI
6636666U60600R060*000006000000000006





AT LOCAL MAX: VELUC = 111. 0 CM/SEC4 HEARING 27. DEGREES


Figure 2.23 FWNO output, f = 15.6 Hz.

























Pi MATRIN FOR TEMPORAL FREOD 17.2 HZ
IN -UU W.*N*. PEAK POWER AT 7.8 HZ
OF 01961399E 0A PICEDATT05EC*CM4*2

37 40 35 3R 35 30 29 31


33 2? 28 22 32 s1 13 IA


3? R4 39 35 27 23 25 33


am IT a 21 16 1 a A II


Is mI 1i 12 ---- I T


20 17 20 mo 22 6 a mI


32 3 31 30 *3 29 31 47


34 25 IO 25 3a 16 14 as


X VELOC


0.250 60.75


.O*e1 91.6P


0.12S 137.50


0.063 275.00


0.0 ******R


-0.063 -275.00


-0.125 -137.50


-0O.16 -91.67


-3 -2 -1 0 I 2 3 4


Oa m i*MMllASAe AMMMI XXM NMMXM NMMEA.




SKOOMBAmeX MXXKXXKKwXsx x ^ XXXam X4 X
MKlX rKMIIMMXNMXMMMXK*KRxMMati )a6
**o^***** ***** a *-**i


M * ** a a Bxx ^*^a f ----g
*L *^* f B* I tl. e.. .. ...-----

** **a*******B=BD***Ba-*- ***---
.** XX* **+*** ****= e ...a---- ..
tp eda it@t wamus twn-- .. ...
.AARAARRAA Ef EXlaAEAXIA.X -----aP




E**R******O.UE=XMMMAXXAAA A XXRE a

SL AL VL 9 CNG DE...


AT LOCAL MAX: VELOC a 91.b7 CMfSEC, HEARING 0. DEGREES


Figure 2.24 FWNO output, f = 17.2 Hz.
























Pw MATRIX FOR TEMPORAL FREDO 10.8 HZ
IN -OB .;R.T. PEAK POIER AT 7*R HZ K VELOC
OF 0.1961J99E 04 PICOTATT*SECSCMN*2

34 36 41 30 31 33 30 29 4 0.250 75.00


33 33 3 o 50 41 20 25 27 3 0.186 100.00


30 42 42 30 28 24 26 33 2 0.125 150.00


35 34 20 31 29 27 26 29 1 0.063 300.00


34 32 25 27 46 27 24 25 0 0.0 *******


33 40 27 27 31 27 27 30 -1 -0.063 -300.00


3J 36 37 30 33 29 20 39 -2 -0.125 -150*00


30 40 46 36 31 26 26 32 -3 -0.168 -100.00



-3 -2 -1 0 1 2 3 4

UUUUUUAMUUUGDaDUUDPNKEKK KNERNESS RKN
BgIg*EMWIB**I***W*I*MWXMKXXMNXXKXXXN
gSEagsE U.UU U K U UU .0aEKES11EAN KEN K
UUUOSEUUSUUUUUUUUUUUDIKAK NEXXSX NNKK
SD ED S U SUE SE SE SEB E KSN EN SE

IIIIIIIIIWlXXIIIXWX IXIHMXXKXMIIXNNM
sIqINIeIIIIIIIIIRNMMMAMMERMMMMMMZW
UsUUEaSSUUSXSSXEm*aNKN:ENN.XXNNEDRNE






ATmi L a MAxx xx x A0 CxS C, loNx M xxxx GME
044I4IIISIEXZIIxx NEXSKX A XXXKNSESNNS








AT LOCAL AM: VELOC 106.07 CH/SEC, HEARING E 45. DEGREES
3D SD K SESSSEN SUE S N KNEES NEESNE NED E


Figure 2.25 FWNO output, f = 18.8 Hz.

























RP MAT0IM FOR TEMPORAL FRED -20.3 HZ
IN -DO *O.0T. PEAK POSEa AT 7.8 HZ K VLLOC
OF 0.1961399E 04 PICO ITT*SECoCMS.2

30 36 34 34 32 40 33 28 4 0.250 81.25


31 33 34 40 37 29 26 27 3 0.180 108*33


38 45 40 34 29 26 28 33 2 0.125 162.50


34 36 27 29 34 29 20 29 1 0.063 325.00


32 36 25 26 34 28 26 27 0 0.0 ***4*00


31 41 30 27 29 29 29 30 -1 -0.063 -325.00


36 35 47 33 33 30 29 34 -2 -0.125 -162.50


39 44 41 35 30 28 29 32 -3 -000.L -100.33



-3 -2 -1 0 1 2 3 4


*05* 00BB000A0005 RM 50N00 I0605 **0 MNN














AT LOCAL MAX V C 32. CMSCC GEARING 160. DEGREES
OOOOUaOOUBUOOMOOOOUIIOOOOS5455SSNSSS
43BABANXIMIMEIMIMIMIIDIMMMEMMMMMMM




lAILAIIIIIIIIIIOIl1XXAMEMM*nKxMN
RMBI0050 5IIII 55 MAISSM5M05500055N5

MIRABAOIWaDMMXXMX*MMNX*MMMMMMMMMXXX
00 00550* N0 A*A*0N00000*00*550ORM*


MMSMIIIIAXXXXNNX XNEMININNNMMMMMARXX
NHAMMBRNaMXXXXXMMNN*INMMMKNMRXNMMXXX
05 OtOOSOS 5* *M ****5055500505500 N*0





AT LOCAL AX: VELUC = 3ZS.DO C0/SCC, OEARING 0 1000 D0G5E5S


Figure 2.26 FWNO output, f = 20.3 Hz.

























PW MATRIX FOR IEMPURAL FRED .21.9 HZ
IN -OU .Re.T. PEFA PDOE34 AT 7.8 HZ
OF 0.1961399C 04 PICO.ATT*SEC*CM*L2

20 36 31 34 35 33 36 20


20 26 7 2-3 2 25 13 12


43 44 35 32 33 JO \30 34




Is1 14 16 Is




16 1I 20 I1 16 1S 5


36 35 39 40 37 33 1 33


24 2I6 32 3 23 2 13 13


0.250


3 0.100


2 0.125


I 0.063


0 0.0


-I -0.063


-2 -0.125


-3 -0.160


-3 -2 -I 0 I 2 3 A

JIIMUSEJNaXINXMmNESN0NgEMXSsISXIXNNX


*XsonMMMaxxxxxxlNxxXXXrXXXRmXXwNXXX


858508.............
MMxxx uxxxMxuxxxxuiu trxx --xues xi
O I I1IXXXXXIIXXIIXxx ..X. KXXX
II IBMM l .IX...I





la.i.. 44.4...i5---a.wmass----
aa:awat*###smanswanasma....


iLiAN0 N"N*""* E s0a0ons5 -.....






LOCAL E 116 II IIIIII BE ING 0. EREE

AT LOCAL M84: VELOC 1363.7 CM/5CC BE*ARING C D. DEGREES


Figure 2.27 FWNO output, f = 21.9 Hz.


VELOC


&7.50


S11.67


176.00


350.00


80*****


-3O50.00


-175.00


-116.67












59











P MATRIX FOR TEMPORAL FkEO t23:4 HZ
IN -OB .R.T. PA PE OER AT 7.e HZ K VELOC
OF 0.1961399E 04 P1C0kAT TSEC5CMEC 2

29 35 32 39 34 31 35 29 4 0.60 93.75


33 33 32 38 41 37 32 30 3 0.l58 125.00


42 39 35 36 45 32 32 38 2 0.125 517.50


30 46 29 30 40 28 26 2B 1 0.063 37B.00


26 44 29 20 35 31 25 23 0 0.0 6**00**


27 35 41 34 38 37 28 20 -I -0.063 -375.00


37 34 3b6 2 6 6 33 34 -2 -0O.25 -187.SD


33 t5 32 34 34 31 42 34 -3 -0.38e -125.00



-3 -2 -1 0 I 2 3 4

NAXNS0MINNSKIISMMSfRRSEXNROEIIEKIRi


IINailnAiXNNMSSSMEUM llSSENlNNASlXNNXl

BAMEANNAmmEEESNmSSESEEMMEMEREUMESSES

ESUDHESESUEMUEEESENxuSSENENNSSESEERE


XRNiOSEllEEMSNIESESELSlliilAI*INEOA
SAEEOEAE MASENESO SESSESE ENANMNE E HDR
gO IEESES NSS AXWAAANSEUWAANAEESEESA
NONIESEESNONEeEENOSSSNOOAAXNNINIIINN

1BAe*max11a11seagme*IIIXr.11:4








AT LOCAL MAX: VELOC N 93.7S CM/SECE BEATING E O. DEGREES


Figure 2.28 FWNO output, f = 23.4 Hz.
























PO MAITI FOR TPMPORAL FEO *2p.0 HZ
IN -UP W..Q.r PEAX PnOER AT 7. Mi
OF 0.1961399E 04 PICUOATT*SEC*CMt*2

32 32 33 49 33 32 32 30


J3 31 32 43 39 55 34 34


36 37 40 o4 36 33 36 39


29 44 33 33 13 31 27 26


26 36 34 30 34 36 26 24


27 31 43 33 33 4* 29 26


36 33 36 36 30 37 36 39


34 55 39 43 32 31 30 34


K VELOC


0.250 I00.00


0.106 I33.33


0.12S 200.00


0.063 00.00


0.0 *******


-0.063 -400.00


-0.125 -200.00


-0.168 -133.33


-3 -2 -3 0 1 2 3 4

xxx3axa3m3 3anagaag3m A itwK6 KK310MM

omaMMWmssxso::agmeaesemagumasa naME:
aImINM5N644.SoIaIanmagagawsamgasaaaa
Iaos an mIasmapemIamsagasagImaxssama
aImasaaaam uma a ammaaaaaaI s aatma ama
smaamaWaaaaaaaaas aaa B aaIaI aaaa as
aalas :namla llnlaaagaaalga 3l a1111a
ROWIIIeIIIIIIIageIIIMERMEMRRMERM
REIIIIIIIIassatWIIIIIIIIRNRNKRRHARAM

mRImmmIMIaIxmsas malw S 9 aII NNxNM. 3 IA
ATLmmCALaAmX.33 10.0a a 3a 0a33 HEIG 0
60 IaaaasIaagnMo6N maMIaIII NNaKaIo
3*050 mmzaanmonMENEmasaaaaaamsWoxxxx

MaxRonnmageagasuageaomagamamnsna
alIIa*aIamsam IImamam*agmaoaoaaIImaXm


aaommoaaxIaaamsasmaagacsmgsaassamgag
atmmammRagImaIIamIMAteggmIagasasagaa



AT LOCAL MAZ: VELOC = 100000 CM/5EC* HEAPING 0. 00EG02EES


Figure 2.29 FWNO output, f = 25.0 Hz.








61

the -3 db contour. Subroutine GREY (Appendix 0) generates the

printed output. The grey levels are produced by overprinting

as shown in Table 3. Linear interpolation is used to determine

the level between given data points.

The first two figures exhibit no power density greater than

-24 db. In Figure 2.15, a peak power density of -2 db is re-

corded. This is obviously due to the first wave listed in

Table 2. The db number array is rounded to 1 db before listing,

but the determination of peak values is based on single precision

floating point arithmetic. Thus the actual local maximum occurs

at 7 = 35.361135* or kx = -0.063 and k = +0.063. The circles of

constant velocity = 30.0 cm/sec are superimposed on the actual

computer output. From the contour plot, it is apparent that the

centroidd" of the 3 db region approximately corresponds to the

actual vector velocity. The superimposed arrows correspond to

the true velocity. The estimated velocity is given below the

grey plot.
Figures 2.16 and 2.17 exhibit only small power levels, but

the twin peaks of Figure 2.18 obviously correspond to the second

and third entries of Table 2. The inner concentric circle repre-

sents the locus of points of 125.0 cm/sec and the outer circle

42.0 cm/sec. This figure illustrates that the wavenumber reso-

lution is about 0.125 cm-1. A comparison of Figure 2.15 with

Figure 2.18 may lead to the erroneous conclusion that the first

wave is 2 db less than the third and fourth. This is a resolu-

tion problem since at 3.1 Hz the power density at the actual

maximum (kx = .087, ky = 0.54) was not computed as it was in the








62
latter case and, instead, we have the peak power density near

the maximum.

In Figure 2.21, the circle corresponds to 75 cm/sec and the

two peaks account for the fourth and fifth entries of Table 2.

The circles of Figures 2.24 and 2.27 both correspond to

100 cm/sec and the peak of these frequencies corresponds to the

last two entries in Table 2. Figure 2.27 also illustrates the

effect of spatial aliasing. The actual wave has a center

wavenumber of k = 0.219, but,due to truncation, the actual

wavenumber is spread beyond the aliasing point of k = 0.25.

As a result the tails of high wavenumbers begin to fold back

onto negative wavenumbers. This accounts for the -11 db level

at the left side of Figure 2.27. The effects of sidelobe leak-

age can be observed in Figure 2.24 above, below and to the left

of the main peak. As these sub-peaks are not much more than

10 db below the main peak, caution must be exercised so that

they are not misinterpreted as additional planewaves.

If the k scale of Figures 2.13 through 2.29 is carefully

observed, it may be noted that there is one more positive term

in each of the k directions than there are negative terms.

This is a characteristic of the FFT which has no counterpart

in one-dimensional frequency series. To understand why this

occurs, consider a band-limited function as shown in Figure

2.30(a). If the time function is sampled at exactly the Nyquist

rate F, the frequency function of Figure 2.30(b) results. If

the transform of the periodic extension of the truncated sampled

function is performed, we have the DFT as shown in Figure 2.30(c).







63


SX(f)



........ ......--------- f

(a)













F
(b)









A ,A

...; --y----------r---------
F
(c)


Figure 2.30 Example of positive and negative frequency terms.








64

If there are N points in the series the first (N/2) + 1 correspond

to positive frequency and the remainder, (N/2) 1, correspond to

negative frequency. Since the DFT of real one-dimensional series

is always symmetrical with respect to positive and negative fre-

quency, the calculated negative frequency may be ignored. In the

wavenumber space, however, this symmetry does not exist and there

is a difference in the number of points on the positive and nega-

tive scales. This apparent deficiency is easily eliminated by

performing a second transform with the spatial directions reversed.

Now the negative frequencies will be computed to (N/2) + terms.

For two spatial dimensions the x and y axes are negated by not

calling NEGATE before transforming to the frequency-wavenumber

domain. If necessary, the x and y axes could easily be negated

separately to obtain the full range in all coordinates.


Summary

This chapter has presented a detailed account of the frequency-

wavenumber space. This should be beneficial in studying the more

complex subjects of multidimensional filtering, Chapter III, and

maximum-likelihood estimation, Chapter IV.

The principal contribution of this chapter is the computa-

tionally efficient direct FFT method of estimating the frequency-

wavenumber power spectral density. This, too, will be used in

succeeding chapters.















CHAPTER III

MULTIDIMENSIONAL DIGITAL FILTERING


This chapter has two purposes. The first is to present a method

of realizing a digital filter (one-dimensional) for an arbitrary

response specification given in the frequency domain. This technique

is needed in the implementation of the maximum-likelihood filter

which is, in turn, a necessary constituent of the maximum-likelihood

estimator. The second purpose is to present a multidimensional

filter which can be used to stop, pass or modify a traveling planewave.

This filter is necessary to reduce the multiwave maximum-likelihood

estimation problem to that of a single wave estimation problem.


One-dimensional Filtering

Recursive filters have an impulse response of infinite duration

and the output of such a filter is a function of both the past output

and the past input. The digital approximation to these is based on

classical analog filter theory and leads to sophisticated design tech-

niques such as Butterworth, Chebyshev, and elliptical filters [133.

Nonrecursive digital filters are approximations to filters which

have a finite impulse response and the output of this type is only a

function of the past input. The principal advantage of the nonrecursive

filter is the much greater degree of design flexibility. Only nonre-

cursive digital filters are considered in this study.

65








66

The digital approximation of nonrecursive filters can be approached

in one of two ways. The first method, the window method, truncates and

samples the impulse response of an ideal filter by using a suitable

truncation window. The actual frequency response is then the convolu-

tion of the ideal response with the transform of the truncation window.

The problem of choosing a suitable window function is the same problem

encountered in forming acceptable power spectral estimates. A draw-

back of this technique is that the optimality criteria must be expressed

in terms of the transform of the window rather than explicitly in terms

of the resultant actual frequency response. Helms [14] shows how the

main-lobe width and sidelobe width can be adjusted in an optimum fashion

using a Dolph-Chebyshev window.

The second approach to nonrecursive digital filtering overcomes

this difficulty for special cases by choosing transition values of the

frequency response so as to minimize the maximum out-of-band response.

This method, called the sampling method, is given in [28]. The syn-

thesis proceeds as follows:

1. the exact frequency response is specified at equidistant points

in the frequency domain for the stopband and passband;

2. the transition points are selected by using computer optimiza-

tion techniques which calculate the maximum out-of-band response as a

function of the transition values. The frequency response is interpo-

lated to any desired degree of resolution by using FFT techniques as

explained later. Computer optimization is used because it is impossible

to obtain general closed form solutions. Unfortunately, it is not

possible to establish general filter design guidelines from [28].








67
Instead one must consult the tabulated results for optimum transition

values. If a filter is required which was not considered in [281, the

costly optimization must be repeated.

The window method and the sampling method are used in subsequent

data analysis; therefore, an example of both for low-pass filtering is

given. In Figure 3.1, a low-pass 16-point filter is given with the

narrowest possible transition width. Only positive frequencies are

shown. The imaginary terms are specified to be zero. The out-of-band

terms are specified to be zero, but appear at -90 db due to plotting

conventions. An inverse FFT is used to transform this function to the

time (impulse) domain. The result is shown in Figure 3.2. The ordering,

characteristic of discrete transforms, should be noted in this figure.

Positive time is in ascending order (9 points), left to right, followed

by negative time (7 points) also in ascending order. The impulse re-

sponse is padded with zeros as shown in Figure 3.3. The interpolated

frequency is then obtained by a forward FFT as shown in Figure 3.4.

For comparison purposes the window (Hanning) method is applied to

the ideal response of Figure 3.1, yielding the specified response of

Figure 3.5. For the specified Hanning window,smoothing convolutionn)

was performed in the frequency domain rather than by multiplication in
the time domain in order to eliminate a redundant transformation. The

impulse response was then obtained as before, yielding the response of

Figure 3.6. After zeros are added and the series is forward transformed

to the frequency domain, the interpolated amplitude response of Figure

3.7 is obtained.













































i
i















s I
o eI r
0 ~ 0


I I I







I I


iI a
















C. uj.u...I..
** (- p
C' o n
f <- c*
c* 0
** i{ I I


I I

1


I I









i i
i I






I I


I I








I I 0




















g


C
v,
w
U,
0

"0
r
to

-)
,
S.-



I-
CL

4.)
I.

2 US
US
'U
'- 0.
t- 3

C r
S5

3
I C
0

US
'U *
*r-
m
c,
S.->

E


c


2
C
S

'-4


C"
s.-
3
01
*-:


























-----C hJ f --* .-- UU -- ,





I
I
I


I !


I 9
I I I




I
I
I I I
I I


I I !
I I


I I

I I I


SI I '
0 0 C r


00
I




1

a










4-
L
I
r
C.





0(-
I




1 : .



1 a0




I S




E









--------------
Sg










70













----,------------------ --------












I I


















- - - - - - - - - - - - - L
j j








I --

5--


I ,U



1 -




;, i n























*---- ------ ---------- ------IA

i I I I

I I S0






41
i I 1

I I,




9 I

I I
I I

* -


I I
I i-





-- ---- ---- -- --- M--- --- - - - -
i I I
I I I 5r
I I I 9 I

I I

















! ..
I iU d i 3
I I l lrl llC
-------- ---. 0
^ = 8 i; 643
B oU,
"323t~%C
ouu-o'U









72








0

i S i S .
S I I I
I I "
a!I I
i i
I S I I a
I I *


I c
I I I I I



I '
5I I I
I i *
I I I
I I I i


-I I I -
I I I






i !
I : L




r
I I I







-r"
I I I


o a. a- a* ( c
* I I I S S Sl w
S S I I S C








0
.I I I I I i
4-i
: I I I I I )
I I I : IC
I (U
5-,
44
I:I :I I:3



I Id
I I
I : I : : La


ClCl II C I 0 Cl C


S I I I I I S S

























I U

* I I


* I I S S I

\I
ISI S 0

I I 0

I I I I 1 0
I I
I


I I~ c

SI I .1 I
I *




I 2 -
r --.-l._ll.- ^. lrrr^ -- l---- __ __ ___ __ -- --- ^ -


i I II
U ,



S *S
I Id








* o

I '


i i
0). 0
,1,
I
s
!__ .-_ __ _- -- - -
S o >> o .> r-0
o n CTt^ o
>^ ^rft n < o
Q O a o o





















0


00
^-^-^-- l~-r~---^ -^ --^l---- --.- ^ ^ -.^ ^ .^
1 ) 1 I t I l I C
I I I I I I I

I I i I i I i I I


SI I ,
i i
I i I i I I S UJ


i Ii I I
I:I I g a.
I II I
l l I I I I I a



I I I
I I I I s I *



I. I I c
I I I I I I I I 4-



Si I I L
I I i
I I I a. !f

I I I I IM I
I I I



Mp4-* ^-Ffc*- FJ~k*!-^- ^-Mr *~ru^ -^-- UBJ *ri--- *: *
S I I I I 0
I I I I


I I 1 I I I s
I 'a C v
I C C
I : 0

I. i : c
I i




i i 1
I 'N *





I I
I 3


I I I I I i

i i11 1 j c>
II I I I I III I
I I I
III II I *
I I la
I lI f r
I I I II I *


-1U-- --- 1---~---- 1MkkJ--^q- ** _ Klc-F----bkh-------11- U U *
P'it o o f o
oP r r c e
SM I j I I i < r I
I I I I I I I I I
* S I I I I I I I










75
















i s 1
II
1 I I 0
I I I I I I! -
sL
i I I I I a)
I i0



| I o L
I I

I a,







__I S S __ I I CS' 0
















aI o o o
O J I I'U Jr














































































ii








----------------- ------ *------ ------- ---


0 a 0 0 0 0 0


I
t




I









77













I I I * I
I I I I
I
I 1 1 1 0
I i I i I SI

II 1
S I i : 3 I
1 -I
i I i
I I I I I .
I S I I lI
I I II
i ** * *1 o
I I I J I
I I "



*I .m
I I I I r

S II I








1 n
I I I I 1fI



S* S II


* SuI I S S I I S
!I




I I o









OIuu aj-








78
The optimum 2-point transition values of the sampling method are

illustrated in Figure 3.8. The impulse response is given in Figure 3.9

and the interpolated response is given in Figure 3.10.

The response of Figure 3.4 which is, in effect, due to a rectan-

gular window is seen to have a sharp transition but a rather large out-

of-band response. The Hanning window method lowers the out-of-band

response, Figure 3.7, at the expense of a wider transition region. The

sampling method provides an additional reduction in the out-of-band

response with no further increase in transition bandwidth. These re-

sults obviously argue in favor of the sampling method, but it must be

emphasized that this technique can only, at present, be applied to a

limited class of filters. The window method, however, can be applied

to any desired frequency response function. The maximum-likelihood

filter of Chapter IV requires a filter which is a function of the power

spectrum of the input data. Since sharp peaks in the spectrum result

in very narrow transition in the filter amplitude response, it was

found that the filter could be successfully approximated by first

limiting the dynamic range of the spectrum and then applying a Hanning

window.

Subroutine AHANN of Appendix D is used to Hanning smooth an arbi-

trary series in the frequency domain. For bandstop or bandpass fil-

tering of a given time series, subroutine BPF may be employed. BPF

automatically smooths the specified transition and returns the filtered

time series. Subroutine SPBPF is an example of a special purpose (i.e.,

fixed transition and fixed number of filter points) bandpass filter








79
which utilizes the sampling method. This filter is used in Chapter V

to pre-condition data for multidimensional analysis.

All of the digital filters developed for this study are specified

in the frequency domain as real functions. In a simulation sense this

is equivalent to zero phase shift. To determine how closely this con-

dition is approximated, all transformations were done using complex

arithmetic and the imaginary interpolated frequency response was

monitored. It was found that for filters of all dimensionality the

imaginary component was consistently at least 30 db lower than the

real component. Consequently, it is reasonable to assume that negli-

gible phase shift has occurred.


Three-dimensional Filtering

In most practical instances sensor data will be available in a

spatially discrete form. It may, as in Chapter V, be temporally dis-

crete also, but, for the moment, this condition will be ignored. Under

these conditions the most general linear filtering that may be applied

to the data consists of separate linear filtering of each channel

(sensor output). This is referred to as generalized multichannel linear

filtering and is illustrated in Figure 3.11. Delay and sum, weighted

delay and sum, and even maximum-likelihood filtering can be put into
this form.

If the spatio-temporal data are not sampled and are treated as con-

tinuous in all dimensions, the most general linear filter is the multi-

dimensional filter. The basic properties of the three-dimensional

filter will be considered. The N-dimensional case is an obvious exten-

sion.







80




s(t,Fr) -----i h s'(t,l)



s(t,r2) -- h(t)2





s(tK) ----- s'(t,rK)





Figure 3.11 Generalized multichannel filter.


=- s'(t,x,y)


Figure 3.12 Generalized multidimensional filter.


s(t,x,y) = h(t,xy)








81
The three-dimensional linear time and space invariant impulse
response h(t,x,y) is the response at time t and position (x,y) to an

impulse applied at the space-time origin. If s(t,x,y) is the input
to such a filter, note Figure 3.12, the output is given as the three-
dimensional convolution of the input with the impulse response,

+0
s'(t,x,y) = /// h(T,c,B) s(t-r,x-a,y-B) drdadS. (3-1)


Furthermore, since s(t,x,y) has a three-dimensional transform

+" -i2in(ft-k x-k y)
S(f,kx,k ) = f;/ s(t,x,y) e dt dx dy, (3-2)


the output can be expressed as


s'(t,x,y) = F~3 S(f,kx,k ) H(f,kx,k )

+0 +i2n(-)
/// S(f,kx,k ) H(f,k ,ky) e df dkx dky, (3-3)


where H(f,k ,k ) is the complex three-dimensional frequency-wavenumber
response of the three-dimensional filter.
The spatial transform of the frequency-wavenumber response yields

i2n(1.r)
H (f,r) = I H(f,k) e dt (3-4)


where r = (x,y) and ( = (kx,k ). The complementary transform is
x y








-2n( ) 82
H(f,t) = / Ho(fr) e di (3-5)


where H (f5,) is an infinite set (with respect to space) of one-
dimensional filters. If we approximate H(f,t) in (3-4) by restricting
Ho to be a sample point operator and limit the number of samples to K,
we have

K -i2 (.(3j)
H(f,) = H (f,5 ) e (3-6)
j=1

The question of obtaining acceptable approximations is readily recog-
nized as a problem in obtaining acceptable estimates of multidimensional
Fourier transforms using truncated data. It has thus been shown that if
spatial stationarity exists, we may transform a multidimensional filter
to K one-dimensional filters. In surveying the literature on the sub-
ject of spatio-temporal filtering, it was found that this technique is
commonly employed in seismology to effect velocity filtering of plane-
waves [29,303.
The filter requirement is usually developed in the multidimensional
domain, transformed to one-dimensional filters, and implemented using
discrete convolution. In this chapter an alternative solution is pre-
sented which consists of digital filtering directly in the multidimen-
sional domain. This method offers substantial savings in computational
effort if the multidimensional FFT is employed. Table 4 compares the
number of operations required for three different, although equivalent,
filters.
As stated earlier, the purpose of the filter, which precedes the

single wave maximum-likelihood estimator, is to pass a single planewave








83
and stop all others. If we assume continuous and infinite data in all
three dimensions and adapt Wiener's least mean square error as a crite-
rion, it can be shown [31] that the multidimensional complex frequency-
wavenumber filter response is given by


H(f,) = Ps(fs)/EPs(f,t)+PN(f,)] (3-7)


where Ps and PN are frequency-wavenumber power spectral density functions
of the signal and noise, respectively. For the required filter the sig-
nal is a single wideband planewave with specified vector velocity


v0 = f/lk2 (3-8)


The noise is the sum of all other planewaves. Thus, the numerator of
(3-7) is zero except when this condition holds and, since the frequency-
wavenumber distribution of signal and noise is unknown a priori, the
maximum uncertainty condition must be assumed, that is, both the signal
and noise are assumed to be flat, and the final filter is

1, t = Vf/lVl2
H(f,l) = {1 (3-9)
0, elsewhere

for all f.
The response pattern is sketched in Figure 3.13. The ray emanating
from the origin is the locus of points I(f,t) = 1 and corresponds to
points for which _0 = ff/|1t2. For all other points in the frequency-
wavenumber space, H(f,) = 0. The resultant filter is a vector velo-
city bandpass filter. It is optimum under the conditions described



























X


-4;




.r-

II



II I

II II


t> t>
*If








85
Table 4
Number of operations required for various


data array sizes


Multichannel
frequency do-
main filtering


13x103

151x103

3x106
4x107


Multidimensional
frequency domain
filtering

9x103

45x103

0.4x106

0.2x107


Multidimensional
discrete
convolution

106

107

108
1010


and its output contains the signal (with no filtering distortion) whose

velocity is vo and any interfering planewave (noise) whose velocity is

V0. In practical terms this means that if two or more waves are with-

in a beamwidth of each other, that is, their vector velocities are

approximately equal, the filter is incapable of separating them and

simply passes the sum as one composite wave. In all other cases, how-

ever, the ideal filter passes only one specified planewave. The remain-

der of this chapter is concerned with the development of digital filters

which approximate the ideal planewave filters.


Implementation
The digital implementation chosen for this study involves approxi-

mating the function of (3-9) by using multidimensional frequency domain

specifications by the window method. Before proceeding, it is important

to realize that we are approximating an idealized response based on


Data
array
size


4x4x32

8x8x32

4x4x1024

8x8x1024








86
continuous data of unlimited extent in all dimensions. As an alter-

native one could reformulate the problem to be one of determining the

optimum linear processor which operates on sampled data of fixed extent

in all dimensions. The relative merit of these two approaches is an

issue aside from the main purpose of maximum-likelihood estimation and

is not investigated here. The frequency domain method was chosen be-

cause it is consistent with techniques previously discussed and appears

efficient from a computational viewpoint.

The principal programming units are divided into four subroutines.

The first three units generate the complex frequency vector wavenumber

response in the frequency domain. The fourth unit utilizes the speci-

fied response and filters a given spatio-temporal series. In addition

to providing the. capability to pass a single planewave based on its

vector velocity, the filters can pass or stop groups of waves based on

their bearing, speed or temporal frequency. This generalization is

provided to facilitate data analysis.


Angle Filtering

The purpose of an angle filter is to pass (or stop) all waves coming

from a particular direction and stop (or pass) all others. The filtering

is to be accomplished independent of the velocity magnitude and band-

width of the planewaves. The relative ease with which this may be car-

ried out stems from the fact that the vector wavenumber t and the vector

velocity have the same direction


v = f4/jil2 .(


(3-10)








87
To specify the frequency response all points lying within the beam

wedge, as shown in Figure 3.14, where o is the specified bearing angle

and is the specified beamwidth, are set to unity, and all other points

are set to zero.

It was demonstrated in the section on one-dimensional filtering

that it is usually advisable to smooth the transition from stop to pass.

This same principle applies to any multidimensional filter. Two-dimen-

sional smoothing need only be applied to each kx-ky plane since the

angle filter response is independent of frequency.


Digital Angle Filter Implementation

Given the values of the complex frequency-wavenumber response only

at discrete points of the ideal response, the complete performance of

the digital angle filter can be determined in one of two ways. The

sampling theorem, as applied to multidimensions, may be invoked to de-

termine the response at every point in the frequency-wavenumber space.

Alternatively, one may numerically interpolate the response utilizing

multidimensional FFT techniques. The latter technique is applied here

as the functions resulting from the mathematical technique do not gener-

ally provide a concise overview of the filtering operation. This is

particularly true when one is interested in examining the amount of

passband ripple and stopband attenuation. This approach is readily rec-

ognized as an extension of one-dimensional, nonrecursive digital filter-

ing design.

Subroutine FA accepts as input the desired bearing (look angle) and

beamwidth and assumes that a 4 x 4 point spatial array, with 2 cm




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs