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
*
* |