ENGN UMR0930'~:Reporft 013-514-1 W.A r) X:: THE USE OF COMPUTER-GENERATED t' PICTURES TO EXTRACT INFORMATION FROM UNDERWATER ACOUSTIC TRANSFER FUNCTION DATA G.; N. lderquist COOLEY ELECTRONICS LABORATORY Department of Eilectrical and Computer Engineering The University of Michigan Ann Arbor, Michigan 481'05 April 1975 Technical Report No. 2-27:. ~ Approved for public release; distribution unlimited. Prepared for OFFICE OF NAVAL RESEARCH Department of the Navy Arlington, Virginia 22217

rvW~Op THE UNIVERSITY OF ML 24A ENGINEERING LIBRARY

SECURITY CLASSIFICATION OF THIS PAGE (When Date Entered) REPORT DOCUMENTATION PAGE BEFORE COMPLETING FORM BEFORE COMPLETING FORM REPORT NUMBER Z. GOVT ACCESSION NO. 3. RECIPIENT'S CATALOG NUMIER 013514-1 4, TITLE (end Subtile) S.. TYPE OF REPORT * PERIOD COVERED THE USE OF COMPUTER-GENERATED PICTURES TO EXTRACT INFORMATION FROM UNDERWATER Technical ACOUSTIC TRANSFER FUNCTION DATA 6. PERFORMING OR". REPORT NUMBER TR227 7. AUTHOR(a) S. CONTRACT OR GRANT NUMBER(e) Gerald N. Cederquist N00014-75-C-0174 *. PERFORMING ORGANIZATION NAME AND ADDRES - 10. PROGRAM ELEMENT. PROJECTX, TAt Cooley Electronics Laboratory AREA 5 WORK UNT NUMBERS University of Michigan Ann Arbor, Michigan 48105 I.I CONTROLLINl OFFICE NAM AND AODDRESS 12. REPORT DATS Office of Naval Research Department of the Navy 1. NUMBeR OF PAGES Arlington, Virginia 22217 212 14. MONITORING AGENCY NAMIE& AoDRESS(I different r Cn.c Itolnd' Off11.) 15. SECURITY CLASS. (of l' r'et).. Unclassified 1 DOIECL.ASSIFICATION/'DOWNGRADING-..CHEDULE IS. DISTRIBUTION STATEMENT r(o. le Rep.). M>4> t <, -.,' Approved for public relea:se-;. distribution unlimited I?. DISTRIBUTION STATEMENT (.1 tie weemdseel'M leok If dif, I e't kei' ReeM) *l. SUPPLEMENTARY NOTES It. KEY WORDS (Centrme n e, eree olk it neoeeM and Dde"lfpy by eblock nmmi,) Underwater acoustic transfer function data Computer graphic-s Underwater acoustic propagation Data display techniques 0. ABSTRACt (CniMue m aen 1ver aid Iof,efee and ider "ed? block m r) A fundamental block to the experimental measurement of the time-varying transfer function of an underwater acoustic channel has been the inability of the scientist to cope with the necessarily detailed and hence large amount of data involved. This research shows that computer-generated pictures are one key to solving this data deluge problem. Data display techniques are developed to scan large amounts of DD,,A,,7 1473 EoITIoN or, NoV r.. oBSOLET /N 010o o014o* 601o I SECURlTY CLASSiFCATION OF TNl PAG Ihr, Da, nered

SECUMiTY CL'A$SSIlICATION OF THIS PAGC (ihM. DOe Conlored) Block 20: multifrequency, bandpass, underwater acoustic transfer function data in the form of computer generated pictures. Transformations are defined to reveal regions of interest as peak.s and edges in the pictures. Thus an observer can follow these regions both in frequency and in time, gaining insights into propagation which are impossible to glean from single-frequency- data. In addition, the techniques permit easy separation of propagation variations due to low signal.-to-noise ratio from variations due to oceanographic effects. These display techniques, applied to transfer function data collected in the Straits of Florida, reveal new insights into underwater acoustic propagation in this locale. 3Cu.TRlY C4A^StCATION 0O THIS PA,.oEI, d.d.. nt.t.-)

ACKNOWLEDGMENTS A research and publication effort of the magnitude represented by this report cannot be solely the effort of one person. An acknowledgment in words to those people who helped in this endeavor is insufficient, but perhaps it will indicate nonetheless my feelings of gratitude. Ted Birdsall served as counselor and pilot for much of the work here. Without his help in the areas of signal processing, underwater acoustic propagation, and communication theory, much of the work could not have been done. Hal Borkin provided the original computer program from which the parallel projection display procedures were adapted and advised on the human perceptual aspects of the computergenerated pictures. Bert Herzog freely gave his opinion on both the computer graphics and my use of the English language. If the reader considers this report well-written, it is due largely to him. (Needless to say, the converse is not true.) Dick Phillips is largely responsible for the human engineering of the display procedures, having made many suggestions to make them more flexible and easy to use. Discussions with L. L. Rauch produced insights on the role of phase as a measure of the travel time in a communication channel. Deep thanks are also due to those people who helped to prepare the work for publication. Carolyn Austin prepared iii

ACKNOWLEDGMENTS (Cont.) the text and Betty Willsher prepared the illustrations, continually amazing me- with their high-quality work while working under tremendous time pressure.. Mary Beth Ogilvie read the first drafts of this report and commented constructively on the presentation. She also confirmed many of the results by repeating the transfer function extraction and display work. Heartfelt appreciation goes to Rosalie Baine whose encouragement and help with the everyday necessities of life has made the period during which this report was written a little more bearable. I wish to acknowledge the support of the Department of the Navy, Office of Naval Research, Contract N00014-67A-0181-0035 (the contract -preceding N00014-75-C-0174) for the work presented here. iv

TABLE OF CONTENTS Page ACKNOWLEDGMENTS iii LIST OF ILLUSTRATIONS LIST OF APPENDICES xii LIST OF SYMBOLS AND ABBREVIATIONS xiii CHAPTER I: OVERVIEW CHAPTER II: PAST WORK AND PLAN OF ATTACK 4 2.1 Review of Underwater Acoustic Transmission 2.2 The Spectrum of Underwater Sound Receptions 2.3 Previous Propagation Research 2.3.1 Theoretical Work 2.3.2 Experimental Work 2.4 The Stochastic, Time-Varying Channel Approach 2.5 The Research Plan CHAPTER III: MEASUREMENT OF THE CHANNEL TRANSFER 23 FUNCTION 3.1 Underlying Assumptions 3.2 The Miami-Bimini Acoustic Test Range 3.3 The Transmitted Signal 3.4 On-site Processing of the Receptions 3.5 Selection of the Data to be Analyzed 3.6 Transfer Function Extraction 3.6.1. Fourier Transformation 3.6.2 Reinsertion of the Carrier 3.6.3 Restriction in Frequency o 3.6.4 Phase Delay Compensation 3.7 Summary of the Properties of S(f, t) CHAPTER IV: DEVELOPMENT OF INFORMATION DISPLAYS 55 4.1 A Brief Introduction to the Display Technique s 1.1.1 Parallel Projection v

TABLE OF CONTENTS (Cont.) Page 4.1.2 Contour Mapping 4.2 Design. Criteria and Constraints 4.3 Converqtions for the Parallel Projection Picturles 4.4j Conventions for the Contour Maps CHAPTER V: ANALYSIS OF THE TRANSFER FUNCTION 75 AM2PLITUDE 5.1 Amplitude Maxima 5.2 Amplitude Minima 5.3 Very Low Amplitudes 5. C-onclus ions CHAPTER VI: ANALYSIS OF THE TRANSFER FUNCTION 89 PHASE 6.1 Reconstructed Phase Angle 6.2 Problems in Correc-ted Phase Due to Noise 6.3 Editing 0 for Minimal Variation 6.4 The Phase Rate 6.5.The Power Spectrum of the Phase Rate 6.6 The Power Spectrum of FrequencyIndependent e Features 6.7 Conclusions CHAPTER VII: INTER-LINE CORRELATIONS 130 CHAPTER VIII > CONCLUSIONS AND RECOMMENDATIONS 136 FOR FURTHER WORK APPENDICES 138 REFE RENCE S 186 DISTRIBUTION LIST 189 vi

LIST OF ILLUSTRATIONS Figure Title Page 2.1 Typical deep-sea velocity profile 6 divided into layers 2.2 Ray diagram of transmission in the deep 9 sound channel for a source on the axis 2.3 A multiphasor model for underwater sound 10 receptions 2.4 A transversal filter model for the under- 19 water sound channel 3.1 The physical layout of the Miami-Bimini 26 range 3.2 The bottom profile of the Miami-Bimini 28 -range 3.3 Sound speed vs depth, Miami to Cat Cay, 28 28-29 November 1961 3.4 Sound ray paths along 25~ 44' on 28-29 30 November 1961 3.5 A CM signal (a) A portion of the modula- 32 ting waveform (b) The resulting CM transmission 3.6 The RMS spectrum of a CM signal for 34 L = 15 and D = 8 3.7 Equipment configuration at Bimini 37 3.8 Complex demodulation of the hydrophone 38 reception 3.9 Block diagram of the signal flow for 40 on-line processing 3.10 Transfer functions of the digital filters 41 used in on-line processing (a) Carrier filter (b) Sideband filter (c) Noise filter vii

LIST OF ILLUSTRATIONS (Cont.) Figure Title Page 3.11 On-line power measurements for 46 25 November 1970 3.12 Measured multipath responses of the 47 Miami-Bimini channel, 25 November 1970 3.13 Deviation in dB between Eq. 3.3 and 54 Eq. 3.4 as a function of frequency 4.1 Example parallel projection of a data 56 surface with 2000 base points 4.2 Example contour map of Fig. 4.1 58 4.3 Five views of the same surface at 65 different rotation and elevation angles. The angles are marked next to each view in the form (rotation, elevation). Because of the half-tone reproduction process used, these views -are more pleasing to the eye than the originals 4.4 Three views of the same data surface 67 showing the effects of different renderings. (a) No crosshatching (b) With frequency markers (c) With crosshatching 4.5 Two views of a relatively smooth data 68 surface, showing the effects of different renderings 4.6 Three views of the same data surface, 70 with increasing vertical variation from left to right 4.7 Example transfer function parallel pro — 71 jection, with coordinate system conventions 4.8 Example transfer function contour map, 73 with coordinate system conventions 5.1 The time-varying power spectrum 76 I s, cf(,t) 1 viii

LIST OF ILLUSTRATIONS (Cont.) Figure Title Page 5.2 Contours of the time-varying power 79 spectrum at 5, 10, and 15 dB down from SMAX 5.3 The relative path loss function 81 5.4 Contours of the relative path loss 82 function at 20, 25, 30, and 35 dB relative path loss 5.5 Contours of the relative path loss 84 function at 25, 30, and 35 dB relative path loss. The 25 dB contour is at the estimated wideband noise level 5.6 Expected and measured cumulative distri- 86 bution functions for those amplitudes below the estimated wideband noise level 6.1 The raw phase j(f, t) with principal 91 value convention 0 to 1 revolution 6.2 The raw phase d(f, t) with principal 92 value convention -.5 to +.5 revolu — tion 6.3 The Riemann Surface of the arctangent 93 function 6.4 The phase e(f, t) as reconstructed 95 by Birdsall's algorithm 6.5 The corrected phase o(f, t) 97 6.6 Contours of the corrected phase 6(f, t) 98 6.7 Path loss and phase angle of S(f, t) 100 during a spot fade on line f = 37 6.8 The track of S(f, t) for f = 36, 37, 101 and 38 in the complex plane for t = 155 to 170 6.9 The probability density function of phase 103 angle as a function of signal-to-noise ratio in dB ix

LIST OF ILLUSTRATIONS (Cont.) Figure Title Page 6.10 90% confidence limits on measured 104 phase as a function of signal-tonoise ratio in dB 6.11 The phase 0(f, t) after manual 106 editing to remove features due to noise 6.12 Contours of the edited phase (f, t) 107 6.13 The magnitude of the phase rate 110 filter transfer function. (a) Linearlinear plot (b) Log-log plot (c) Freauencv-period conversion nomogram 6.14 The phase rate 0 (f, t) of the cor- 111rected 0 (before editing). (a) At 112 200 rotation (b) At 0~ rotation 6.15 The phase rate 0(f, t) of the edited 1140. (a) At 200 rotation (b) At 0~ 115 rotation 6.16 The phase rate obtained by setting 117 e(f, t) from Fig. 6.14 to zero wherever the signal-to-noise ratio fell below +3 dB 6.17 Two-dimensional power spectra of the 121 phase rate 8 (f, t) from Fig. 6.15 6.18 Coherent power spectrum of the phase 123 rate along the time axis 6.19 Incoherent power spectrum of the 125 phase rate along the time axis 6.20 Five point moving average of the 127 coherent power spectrum of Fig. 6.18 7.1 The sliding correlation of S(31, t) 132 with S(f, t) 7.2 The sliding correlation of the ampli- 134 tudes of S(31, t) and S(f, t) x

LIST OF ILLUSTRATIONS (Cont.) Figure Title Page 7.3 The sliding correlation of the phase 135 portion of S(31, t) and S(f, t) xi

LIST OF APPENDICES Page APPENDIX A: THE SOFTWARE ENGINEERING OF THE 138 DISPLAY PROCEDURES APPENDIX B:- DISPLAY PROCEDURES FOR PARALLEL 149 PROJECTION VIEWS OF DATA SURFACES APPENDIX C: CONTOUR MAPPING PROCEDURES FOR: 159 DATA SURFACES APPENDIX D: BIRDSALL'S PHASE RECONSTRUCTION 175 ALGORITHM APPENDIX E: THE PROBABILITY DENSITY ON PHASE 178 AS A FUNCTION OF SIGNAL TO NOISE RATIO APPENDIX F: THE MAGNITUDE OF THE PHASE RATE 182 FILTER TRANSFER FUNCTION xii

LIST OF SYMBOLS AND ABBREVIATIONS a. The received amplitude from sound ray path i B Kailath's channel bandwidth spread extent, or the intercept of the channel phase delay compensation line c The speed of sound CM Complementary-phase modulation C-power The carrier power measured on-line at Bimini c(t) The complex-valued time series which is the output of the carrier extraction filter d Duration in seconds of a pseudorandom sequence digit D The number of cycles of carrier per pseudorandom sequence digit dB Decibels di The length of sound ray path i E(x) The probabilistic expectation of x f Continuous frequency, or the frequency index for sampled data Af The spacing between adjacent lines in the channel probe signal fc The carrier frequency of the channel probe signal FFT The fast Fourier transform algorithm fp(x) The probability density function on the power of Gaussian noise H(f) The transfer function of the phase rate filter xiii

LIST OF SYMBOLS AND ABBREVIATIONS (Cont.) H(f, X) The Fourier transform of the channel impulse response h(t,X) Hc (f) The transfer function of the carrier extraction filter HN(f) The transfer function of the noise extraction filter H (f) The transfer function of the sideband extraction filter h(t, X) The channel impulse response Hz Frequency in Hertz Im Imaginary part, or imaginary part of j The square root of -1 L Kailath's channel time spread extent, or the number of digits in the pseudorandom modulating sequence L(f, t) The relative path loss at frequency f and time t M The slope of the channel phase delay compensation line N-power The noise power measured on-line at Bimini Re Real part, or real part of rev Revolutions as an angle measure r(t) The representative time response of the channel to the channel probe signal s(t) The complex-valued time series which is the output of the sideband extraction filter S(f, t) The sampled transfer function of the channel S(k, t) The discrete Fourier transform of s(t) xiv

LIST OF SYMBOLS AND ABBREVIATIONS (Cont.) S(N, c) The spectrum of observable underwater sound-receptions under sound speed structure c SMAX The maximum amplitude in S(f, t) S-power The sideband power measured on-line at Bimini t: Continuous time, or the time index for sampled data T The period in seconds of the channel probe signal u, u' The transform of frequency f in the twodimensional discrete Fourier transform of v Signal-to-noise ratio VI V', Va The transform of time t in the two-dimensional discrete Fourier transform of 6 A Used as an operator to mean "change in," or the time between successive data sets in the on-line processing at Bimini, or the time difference between successive multipath arrivals 0, 9(f, t) The phase of S(f, t) ranging over the real line 6. The phase change in sound ray path i due to reflections at boundaries e(u, v) The coherently averaged, two-dimensional discrete Fourier transform of the phase rate e' (u', v') The incoherently averaged, two-dimensional discrete Fourier transform of the phase rate ea (0, va) The five point moving average of the coherently averaged power spectrum 8(0, v) p The carrier power scaling factor used in computation of S(f, t) xv

LIST OF SYMBOLS AND ABBREVIATIONS (Cont.) p2 The sliding correlation of S(31, t) with S(f, t) 2 The sliding correlation of the amplitude of a S(31, t) with the amplitude of S(f, t) p2 The sliding correlation of the phase of p S(31, t) with the phase of S(f, t) iT The propagation delay time for sound ray path i X, P(f, t) The measured phase of S(f, t) in a one revolution range w Angular frequency <<x, y>> The sliding correlation of time series x with time series y xvi

CHAPTER I OVERVIEW In the past, the propagation of underwater sound has been studied mainly in terms of sound rays or mode structures derived from deterministic physical considerations. Once salinity, temperature, density, and boundary profiles are known, underwater sound propagation between two points can be calculated, provided sufficient computer time is available. The difficulty with this approach is that many of the required physical parameters are difficult or impossible to measure with sufficient accuracy. An alternate approach is to treat the ocean as a probabilistic, time-varying information channel. Underwater sound propagation between two points may then be characterized by the probabilistic, time-varying transfer function of this channel. Past work with single frequency transmissions hints that a small number of distinct channel models, each with slowly varying parameters, should be sufficient to characterize the transfer function of a large number of interesting underwater sound paths. With the advent of small and inexpensive computers, it becomes economically possible to dedicate a machine to collect and characterize transfer function data over long periods of time. Thus it is now feasible to use automation to compile a catalog of transfer functions and 1'

2 associated probabilities for many underwater sound paths. Before any such automatic cataloging effort is undertaken, exploratory measurements of transfer functions over a wide passband (as opposed to a single frequency) will be needed to determine appropriate types of channel models. A major hindrance to this exploratory research has been a lack of tools to handle the deluge of data involved. Past explorations have been confined to one or two frequencies to reduce the data analysis problem to manageable proportions. The research reported here has shown, however, that the application of computer graphics to analyze large amounts of multi-frequency data yields significant insights into the wide passband underwater acoustic channel. In an exploratory experiment, a periodic probe signal was transmitted underwater between two sites 78 km apart. The reception was processed to yield 16,000 samples of the channel transfer function over a 50 Hz bandwidth and a seven-hour period of time. Efficient, interactive computergraphical display procedures were developed and employed as data handling tools for these samples. Pictures were drawn of statistics of the amplitude and phase, and of the inter-line correlation of the transfer function. Wide passband (or wide bandpass as an adjective) is used here both to mean somewhat wider bandwidth than that of commonly-used acoustic transducers, i.e., from onetwelfth to one-half octave, and to emphasize that this study deals with bandpass as opposed to lowpass signals.

Analysis of the resulting pictures (1) indicates a number of features which prospective channel models must include, (2) has already caused a reevaluation of the parameters previously used to characterize single-frequency receptions, and (3) suggests additional methods for processing transfer function data in further experiments.

CHAPTER II PAST WORK AND PLAN OF ATTACK 2.1 Review of Underwater Acoustic Transmission Only those aspects of underwater acoustic propagation germane to this study are presented here. A more detailed treatment may be found in the excellent book by Urick [1]. The propagation of sound in a marine environment is an appreciably more complex phenomenon than information transmission through more commonly studied channels such as telephone cables. The propagation medium is dispersive, inhomogeneous, nonisotropic and time varying. Moreover, the boundaries of the medium are rarely well known and vary with time. For these reasons, underwater sound propagation remains a topic for research. The speed of sound in water is a function of temperature and pressure. Since temperature and pressure vary with depth, the speed of sound in the ocean is a function of depth. Descending from the surface, the speed of sound drops with increasing depth (due to decreasing temperature) until a minimum is reached. With further increases in depth, the temperature remains constant and the increasing pressure causes the sound speed to rise again. The plot of sound speed versus depth is called a sound speed profile. An exemplary sound speed profile 4

5 (after Urick [1]) is shown in Figure 2.1. With such a profile, sound in the thermocline layers (Fig. 2.1) tends to refract downward toward regions of lower sound speed. Conversely, sound in the deep isothermal layer tends to refract upward toward regions of lower sound speed. Measured sound speed profiles do not behave as well as in Fig. 2.1. The thermoclines are usually irregular due to temperature perturbations near the ocean surface. Occasionally, these perturbations change the sound speed profile sufficiently to cause the phenomenon of ducting, where sound is carried for long distances without appreciable energy loss and with no interaction with the ocean bottom. Ducting occurs anywhere in the thermocline layers where the sound speed profile reaches a local minimum. Sound exhibits normal wavelike behavior when it encounters a change in the density of the propagation medium. If an underwater sound wave encounters the ocean surface, it is almost completely reflected with only a minor energy loss due to scattering; it also experiences a 180-degree phase change. The situation at the ocean bottom is not so simple as at the surface. The bottom ranges from thick mud to hard packed sand or rock, each with different layering, each with different density. Experiments show sound waves to be reflected and scattered as well as absorbed when they encounter the bottom. Moreover, due to

6 Surface Velocity of sound, fps layer 4,850 4,900 4,950 50 S easonal hermoc 1 ine Main thermocline 3,000 Deep isothermal layer 6,000 9,000 12,000 Fig. 2.1. Typical deep-sea velocity profile divided into layers

7 the inhomogeneity of the bottom material, no simple statement can be made about the phase changes in a wave when it encounters the bottom. Sound propagating underwater is thus affected by refraction in the medium and by reflection and scattering at the boundaries. In a shallow region with no isothermal layer, sound may propagate between two points by bouncing time and time again off the bottom, "hopscotching" its way across the ocean floor. Many of the physical parameters which govern underwater sound propagation vary with time. The ocean surface varies rapidly with time, and sound reflected from the surface abquires a rapidly time-varying character. Temperature and density vary more slowly. Small perturbations in temperature, density, and salinity have been observed to propagate as waves in the thermocline layers. These perturbations, called internal waves, change the speed of sound at most 0.01 percent. Internal waves have periods ranging upward from five minutes. Diurnal heating and cooling affects sound speed in the surface layer. Tides and seasonal variations in currents and prevailing winds affect propagation on a time scale from days to years. Time-varying biological and man-made noise sources also cause changing signal-to-noise ratios in underwater sound receptions. Physicists have studied underwater acoustic propagation

using two main approaches: mode theory and ray acoustics. Mode -theory views the ocean as a waveguide, using partial d.ifferential.equations to describe the sound intensity.in the ocean as a function of its boundaries. Mode theory has proven useful mainly for;shallow water and small enclosed volumes (-such -as harbors). Conversely, the ray acoustics model of underwater sound propagation largely ignores interactions with boundaries. It has proven useful for longer distances and deeper water. This is the situation of the exploratory expe-riment used for this study (Chapter III). In ray acoustics, the locus of points formed by the normal to-the.sound wave front as it propagates is called a sound ray path. Examples of sound ray paths are shown in Figure 2..2, where each ray path is indexed by its angle of launch from the sound source. Between two geographical points, a ray propagated at a larger angle from the horizontal travels a longer distance, but at higher speed, than a ray propagated at a small angle from the horizontal. As shown in Fig. 2.2, the sound rays present at -a receiving hydrophone may have been propagated from the source at several different launch angles.- The propagation delay., and:thus the arrival instant, for each ray is a function of the distance traveled and the sound speeds encountered along its path. Propagation by more than one path is termed:multipath. A single pulse of transmitted

Range in Miles Sound Velocity,fps Ln C) U1 C> C)~~~~c C~ Hn O H 12. 2 0 0n o ulo0 12.200 110 70 2,000 1 1 05.19 in ~ ~~~~~~ I- *II*'I I 04 a! 12,200~~~~~~~~~~5.1 Fig. 2.2. Ray diagram of transmission in-the deep sound channel for a source on the axis

1') sound may arrive at the receiver as a time spread sequence of distorted, weakened pulses if propagated through the ocean via multiple paths. Frequency spreading also occurs in underwater sound propagation. Interaction with the surface, interaction with internal waves, and nonlinearities in the medium modulate the propagating sound wave in both frequency and phase. A sound wave initiated at a single frequency arrives at the receiver as a band of frequencies. Typically, the bandwidth of the modulating processes is narrow so that the modulated sound wave may be modeled as a multicomponent phasor as shown in Figure 2.3. A strong phasor - modulation process I., phasors transmitted nhasor Fig. 2.3. A multiphasor model for underwater sound receptions

11 at the transmitted frequency combines with weaker phasors of time-varying amplitude and phase. Since multipath propagation occurs, underwater sound receptions may be modeled as the sum of many phasors of the type shown in Fig. 2.3. The resultant amplitude and phase of the reception are susceptible to cancellation effects between sound waves which traveled via different paths. Small changes in the physical parameters governing the individual propagation delays can thus lead to large changes in the resulting observable amplitude and phase. 2.2 The Spectrum of Underwater Sound Receptions The spectrum of, underwater sound receptions at time t' = t + ANet is very nearly the spectrum at time t translated in frequency. To see this, consider sound propagating between source and receiver over a number of different paths. The sound speed c in the medium is a function of position as well as time: c = c(x, y, z, t) (2.1) For each distinct path there is a received amplitude ai, a propagation delay time T., and a phase change 1 1 ei due to reflections at boundaries. These three parameters are all functions of the sound speed c. At time t,: c has a particular spatial distribution which

12 all] ows certalin pa:t]ths to sultpport: sound tprop)acation between source and receiver. Let P be the set of these propagation paths Pi P = {Pi' i = 1, 2,..., N} for some N (2.2) Let the time delay T. be 1 d. Ti = (2.3) Ci where d. is the length of path Pi and ci is the average sound speed over Pi. The spectrum of the reception from path i, Si(:, c), is a function of the reception parameters ai, - and. ]. +je Si(w, c) = ai e (2.4) and the spectrum of the observable reception, S(., c) is by superposition and Eq. 2.3 d. N j +j Ci S(,c) = L a. e (2.5) i=l

13 At some later time t', the speed of sound will have changed to c' = c + Ac (2.6) The medium will now support a different set, Q, of propagation paths between source and receiver. If Ac is sufficiently small and if there are no paths which changed character (e.g., intercept the ocean surface under c but do not: intercept the ocean surface under c' ), then Q will contain N paths. Further, these paths will be geographically close to the paths in P. Under these conditions', the amplitude and boundary phase change parameters will be nearly unmodified: a. = a. (2.7) 0' = 0. 1 1 To obtain an expression for T!, assume that d has not changed appreciably for path i: di = di (2.8)

14 Now di = Ci. * T (2.9) and d1 C1 T1 (2.10) so that by Eq. 2.8, C. T. = c! (2.11) Ci i (2.11) Assuming that Ac is the same for all paths, solve Eq. 2.11 for T! using a series expansion: iT = T. ~ = 7.1 - + higher order terms 1i 1 $~:C C. (2.12) The coefficient?i on the first order term would be one if Eq. 2.8 truly holds. The fact that Eq. 2.8 is only approximate can be compensated by the choice of a. Substituting Eq. 2.7 and Eq. 2.12 into Eq. 2.4 yields 1 c + jei CC Si(!i9, c') = ai e (2.13)

15 Grouping in the exponent yields i ~ i 1 - + j- A S.(a), c') = a. e (2.14) i ( ( acl ) c ) Thus the spectrum of a path under c' is a frequency translation of the same path spectrum under c. If all the path sound speeds ci are nearly the same, then by Eq. 2.5, S(w,c') = S ( (~1 -: c ) c (2.15) 2.3 Previous Propagation Research 2.3.1 Theoretical Work. The propagation of underwater sound is governed by a set of nonlinear partial differential equations. They may be solved numerically to determine the underwater sound propagation between two points. First, one computes the travel time, phase delay, and power density for each sound wave to reach the specified receiving hydrophone. These are then summed via Eq. 2.5 to compute the observable amplitude and phase as a function of frequency. There are some inherent

16 pitfalls in this approach. The solutions are only as good as the specified boundary values for the differential equations. These boundary values are the surface profile and the bottom profile and composition between the source and the receiver. Another problem is knowing the temperature and density as a function of depth and horizontal range. These are nearly impossible to measure accurately and with sufficiently fine resolution. These problems arise in both ray path and mode structure models. Even if the foregoing problems were resolved, large amounts of computer time are required for all but the shortest paths. Thus: results are usually outdated even before they are obtained as the ocean will have changed during the computation period. 2.3.2 Experimental Work. Experiments have been conducted to measure underwater sound propagation directly, mostly by measuring the ocean's impulse response. Usually the impulse is generated by an exploding charge or through the use of an electric spark discharge. However, these impulses do not provide the repeatability in timing and emitted frequencies so necessary for well-controlled measurements. Moreover, the high pressures involved invoke nonlinearities in the medium. Moving platforms are commonly used for underwater sound propagation measurements. Underwater sound wavelengths are quite short, e.g., 4.4 meters at about 400 Hz.

17 Thus even small platform movements can introduce unacceptable measurement errors. Moving platforms also lead to doppler shifting of the frequency of the received transmissions, introducing another source of error. Many previous propagation measurements have been of a survey nature. Researchers have made a few measurements of a small number of parameters over short time periods at geographically separated points. These data, while useful, are not the long term, highly controlled measurements needed to investigate the underwater sound transfer function in detail. 2.4 The Stochastic, Time-varying Channel Approach Consider the underwater sound propagation path as a probabilistic, time-varying information transmission channel. The emphasis then shifts from the physics of underwater sound to the study of the propagation channel using communication theory. The advantages of this approach are that: Extensive and expensive measurement of all of the physical parameters of the path is not required. One can maintain excellent control over a few selected experimental variables (e.g., frequency and timing of channel probe signals). One can relate observed propagation characteristics to pastobservations.

18 If the ocean surface modulation action is set aside and the sound pressures are small, one approach is to characterize the underwater channel as a filter. Previous work with single frequency transmissions suggests that a small number of different filter types, each with slowly perturbed parameters, may prove sufficient to model all types of propagation observed at a particular site. For example, each- different filter could correspond to a particular form of the sound speed profile, with probabilistic parameters expressing deviations from the norm. The transfer function of the filter as a function of frequency and time then gives the input/output properties of the underwater acoustic channel. One specific channel filter model is the tapped delay line, or transversal filter. This realization is quite intuitive because the delay between taps corresponds to the time delays between successive arrivals at the receiver. Such a filter is pictured in Figure 2.4. The tap gains gi are complex-valued so as to model the unknown phase changesoccurring at reflecting boundaries. The tap gains are also random processes which change slowly with time to model the observed slow changes in propagation. The length of the delay line corresponds to the number of significant arrivals. Changes in the number of arrivals indicate changes in propagation mode (i.e., in the sound speed profile) and require the use of a physically different

19 Complex Gains Added Noise Delays Superposition'ransmission Noise 2r 2 0 +W-+ Reception. n nn Fig. 2.4. A transversal filter nodel for the underwater sound channel

20 filter, although the conceptual filter model would remain the same. Even though it is intuitively appealing, the transversal filter may not be a good model of underwater sound propagation. Other models may be better. All proposed models should be assessed by testing their capability to explain observed propagation. An iterative approach to channel model selection is impeded by the very large amount of data to be handled, both in evaluating existing channel models and in deducing new channel models. Evaluation is the easier of the two. It involves computing a normalized error measure and comparing it with previous results. All computations are readily performed by the digital computer. It is in deducing new channel models that the most substantial problems occur. To guarantee sampling above the Nyquist rate, as many as twenty million data points per day are collected in typical wide passband underwater sound propagation experiments. Preprocessing this data to reduce its bulk for human consideration requires that the preprocessor employ a paradigm of the information expected in the data. Such preconceptions can unwittingly lead to discarding important information which does not match the paradigm. Thus it is desirable first to present the raw data for the researcher's inspection. After he has gained confidence that his ideas of information are reflected in the data, he can more wisely guide the reduction of bulk.

21 In general, information in raw underwater sound propagation data takes two forms: repeated patterns of behavior, and large-scale departures from previous behavior (events). To discern patterns and events in so many raw data if they are in tabular form is a formidable task. The human eye, however, is facile at discerning patterns and events in visual presentations, particularly if they can be viewed from different angles. Hence computer graphics has been successfully employed in this study to recognize patterns and events in raw underwater sound propagation data. 2.5 The Research Plan At the onset of this study, high quality experimental data in the form of time series were available from a previous wide passband underwater sound propagation experment [6]. A review of that experiment is given in the first part of Chapter III. To use these data for this study meant that they must be Fourier transformed. Efficient, human engineered computer graphical display procedures would have to be developed and verified and then applied to the transformed data. With these points in mind, the following research plan was developed: (1) Process the wide passband receptions from [6] to yield statistically reliable samples of the transfer function over a selected frequency band.

22 (2) Investigate and develop several raw data display procedures. (3) Display pictures of the transfer function surface and various transformations of the transfer function surface to determine the suitability of the graphical method, and to recognize and extract features in the data.

CHAPTER III MEASUREMENT OF THE CHANNEL TRANSFER FUNCTION The underwater sound propagation data to be analyzed were collected during November 1970 at an acoustic test range spanning the Straits of Florida. Analyses of the signals and receivers may be found in [8], [9], and [4]. Much of the information concerning the test range is drawn from [10]. 3.1 Underlying Assumptions Several assumptions underlie the data collection experiment. In general, they concern the effects of the ocean on the received signal. At sufficiently high sound intensities, the nonlinearities in the differential equations governing underwater sound propagation play a significant role. By using low-transmitted power levels, these effects can be ignored. When using low-transmitted power levels, the receptions will be, in general, below the noise level measured in a 1 Hz band at the receiver. Consequently, long time-constant (i.e., narrow-band) filters must be used to extract the signals of interest. However, the time constants of these filters must not be so long that the ocean transfer function will change significantly during the time the filter develops its output. Previous measurements [10] have indicated that propagation changes 23

24 only slightly over a period of a few minutes. Thus filter time constants of 100 to 200 seconds are utilized. With such time constants, the receiver can improve the signal-to-noise ratio of the receptions by approximately 20 dB relative to a 1 Hz band before detailed analysis of measurement error introduced by variations during measurement becomes necessary. During certain portions of the year, the acoustic paths through the test range appear to interact significantly with the surface. In the Straits of Florida the upper limit on the periods of those surface waves having significant energy is approximately 10 seconds. Thus surface modulation sidebands will appear at least 0.1 Hz away from a transmitted signal line. The bandwidth of any signal line extraction filters used in the receiver should then be less than 0.2 Hz. Investigation has shown it is possible to obtain an "instantaneous" measurement of the transfer function of the test range. Kailath [2] investigated the problem of making measurements on doubly spread (i.e., spread in both frequency and time) channels. In his theory, the severity of channel spreading is measured by the product B.L. B and L are defined in terms of the channel's impulse response h(t,X) and its Fourier transform H(f, ):

25 B = max Bandwidth of h(t,k) mx k Smallest fo0H(f,X) = 0, f > f0 (3.1) L ='max Duration of h(t,k) t max {Smallest Xo1h(t,X) = 0, X > (3.2) If B'L > 1, the channel is called overspread and only ambiguous average measurements may be made. If B*L < 1 the channel is called underspread and instantaneous measurements of h(t,X) may be made. Veenkant [3] showed by actual measurement that B-L = 0.36 for the test range. Hence instantaneous measurements of the test range can be made. 3.2 The Miami-Bimini Acoustic Test Range The Miami-Bimini range, illustrated in Figure 3.1, is a facility of the Acoustics Group of the Rosenstiel School of Marine and Atmospheric Sciences (RSMAS) of the University of Miami. It extends across the Straits of Florida from Miami, Florida, to Bimini, Bahamas.

26 25e Depth in folhormn 50' ~ |,7i MARINE LAB 20o, BIMI I VIRGINIA K-EY 1 I1 6 4N II IS.BIMI I " 10' 280 3 20'4 30' Id "' A T L A N T I C, Fig. 9.1. The. Byic'1421~1,, 10' 8O" 520' 40 30 20' 10' Bimini range The transmitting site is located at Fowey Rocks (Point 2, Fig. 3.1), approximately 19 km from the RSMAS laboratory, and is connected to the laboratory by telephone lines (Point 1, Fig. 3.1). At Fowey Rocks, a bottom mounted acoustic transducer is situated in 22 m of water at the focal point of a 7.3-m compliant tube parabolic reflector. The transducer has a maximum source level of 120 dB/ubar at one yard with a nominal bandwidth of 100 Hz. The source level for the data acquisition parab~oli relco/h rndcrhsamxmmsuc I nu'lr —,nri/hir;f-nn=vadwt anoia bad dt

27 experiment was 110 dB/ubar at one yard. The 30-degree beamwidth is directed toward Bimini, a distance of 78.6 km (43 nautical miles) across the Straits. Two receiving sites were used in the experiment. At the first site, a bottom-mounted hydrophone is located in 300 m of water (Point 4, Fig. 3.1), approximately 12 km (7.5 nautical miles) from the source. The reception from this hydrophone is transmitted by cable to the RSMAS laboratory. The second receiving site is located off Bimini (Point 3, Fig. 3.1). There, two bottom-mounted hydrophones are cable connected to the Lerner Marine Laboratory on Bimini. These hydrophones are located at 30-m and 366-m depths, 1.8 km and 3.6 km offshore, and are referred to as the shallow and deep hydrophones, respectively. The bottom profile of the Miami-Bimini range is illustrated in Figure 3.2. A shelf extends out from Miami about 28 km to a depth of 400 meters, followed by a sharp drop-off to a depth of 800 meters. Fifty kilometers beyond, the Grand Bahamas Bank rises abruptly from this depth. Sound speed profiles for the range, obtained in November 1961, appear in Fig. 3.3. The profiles show This is approximately 300 watts acoustic, concentrated in a beam covering only 4% of the solid angle of a sphere.

28 60.00' 79'50' 79 40?'30 7?9"20' w. 0 LONG 200,:' 400 z k 600 0 Bo800 (VERTICAL SCALE 37 HORIZONTAL) 0 O10 20 30 40 N. MILES a~~~ -___. I (ONE TO ONE SCALE) Fig. 3.2. The bottom profile of the MiamiBimini range NAUTICAL MILES 0 3 20 30 40 - a I i I I I I i I I 5 I I WEST LONG.' DO00 7-.7790'' 7."' z'7J0' _ V. 3. 7..2. 0 o ~ Y, o 0, oc' ~ o~~ P o~~~~~~~'STATION POSITION SOUND SPEED SCALE SEC. L j, o~...'l... n- - "~"0'/E~'50 " "2 1~"... $00 Fig. 3.3. Sound speed vs depth, Miami to Cat Cay, 28-29 November 1961

29 show a constant velocity surface layer which extends to about 100 meters in depth, underlaid by a thermocline of negative velocity gradient. The velocity gradient becomes increasingly negative as the Florida shore is approached. Figure 3.4 shows a sound ray diagram corresponding to the sound speed profiles of Fig. 3.3. This diagram shows that depending upon the ray launch angle, sound propagates from Miami to Bimini by reflection from the surface and bottom, by refraction and reflection from the bottom, and by refraction and reflection from the surface. The refraction/surface-reflection mode of propagation does not appear in ray diagrams computed from sound speed profiles obtained during the spring and summer. Although not indicated by Fig. 3.4, occasional subsurface ducting has been suspected during previous propagation measurements. 3.3 The Transmitted Signal The transmitted probe signal consisted of a carrier wave at 420 Hz modulated by a linear maximal pseudorandom binary sequence. The phase of the carrier was shifted either +45 degrees or -45 degrees depending upon the value of each binary digit in the modulating sequence. Such a signal is termed a CM signal [4]-. A portion of such a

-20 400- I 20*~~~~~~~~~~~0 Fig. 3 4. Sound ray paths along 250 44' on 28-29 Nov. 1961

31 CM signal is shown in Figure 3.5, where f = carrier frequency in Hertz c d = duration of the sequence digit in seconds D = number of cycles of carrier per sequence digit. The parameter D is chosen to be integer valued. Thus the CM signal is periodic with period LD T = Ld = D f where L is the number of digits in one period of the modulating sequence. The CM signal has a number of advantages for'use as a channel probe signal: It has essentially a sin(x)/x line spectrum where the inter-line spacing is controlled by the parameter T, and the location of the spectral zeros is controlled by the parameter D * It has a small peak-to-average power ratio, a desirable consideration with peak-power-limited underwater acoustic transducers. * It is easily generated digitally, resulting in highly reliable timing and purity of waveshape. * Approximately half the total transmitted power is at the carrier frequency, yielding a signal useful for both single frequency and broadband measurement

+1 L dD/fo (a) 1j 1/fe (b) Fig. 3.5. A CM signal (a) a portion of the modulating waveform (b) the resulting CM transmission

33 The RMS power spectrum of a typical periodic CM signal is shown in Figure 3.6. In addition to the advantages cited above, the wellknown autocorrelation properties of linear-maximal pseudorandom sequences (cf. [4]) allow the received CM signal to be correlated with a stored reference to produce the so-called multipath response of the underwater channel. The multipath response is a plot of received signal versus time; it graphically exhibits the multiple signal arrivals due to path diversity in the ocean. The multipath response corresponds approximately to the impulse response of the channel except that the measurement is made with a periodic band-limited signal. Measured multipath responses of the test range will be presented later in this chapter. For additional information on multipath measurements using CM signals, consult Veenkant [3]. The specifications of the probe signal used in the exploratory experiment are: f = 420 Hertz D = 8 cycles/digit L = 63 digits with a resulting digit duration and period of: d =0.019 seconds T = 1.2 seconds (exactly)

a\ H \ r~~~~~~~~~~~~~d~~T o ~~~01 U cl/ ii ~ric) v) o L(I E- - + C <l, c, Uo a__ _ _ r_ II -- -- "O * ~T.......................W ~~~~~~~~~~~~o. ~~~~~~~~~~~~~k (d Ad o Q, k I'~......' —..... (D a~ ~~~~~....I~r-II _ -....C.....r- -[...... _.I - /m O s c. _H m

35 The spacing Af between adjacent lines in the transmitted spectrum is Af = 1/T = 5/6 Hertz The bandwidth of the main lobe of the transmitted spectrum (between the first spectral zeros adjacent to the carrier) was chosen to be 105 Hertz to match the response of the transducer. There are 127 spectral lines lying within the main lobe transmission bandwidth of 105 Hz. Because the multipath response measurement is made with a periodic signal, it is liable to aliasing errors when the duration of the multipath response is more than T seconds. To lessen the possibility of such errors, T should be chosen sufficiently long that all significant signal arrivals fall within one-half of the signal period. Extending T to 1.2 seconds allows some safety factor, as the multipath spread of the Miami-Bimini channel is observed to stay below 600 milliseconds except on rare occasions. 3.4 On-site Processing of the Receptions During the November 1970 experiment, preliminary online data analyses of the receptions were performed. The quantities measured on-line were: signal power, noise power, and carrier angle,

36 * forward scattered surface reverberation, and ~ multipath response. These measurements were obtained using digital signal processing techniques implemented in a small digital computer located at each receiving site. Each of the measurements was computed and recorded every 99.6 seconds on the basis of 64 1.2-second periods = 76.8 seconds of reception. (The remaining 22.8 seconds were used solely for computations.) After each 76.8 seconds of reception was digitized, accumulated, and filtered, various data were recorded on digital magnetic tape for later off-line processing. The equipment at Bimini is shown in Figure 3.7. A frequency standard stable to one part in 1010 drives a frequency synthesizer to produce a 1680 Hz (equal to four times carrier frequency) clock signal. The 1680 Hz signal controls the analog-to-digital conversion circuitry on the computer. The computer demodulates the reception sampled by the A/D converter from a real-valued, bandpass signal centered about 420 Hz to a complex-valued, lowpass signal centered about dc (i.e., about zero Hz). The digital processing mirrors the conventional analog in-phase/ quadrature processing shown in Figure 3.8. No information is lost in this "complex demodulation" process, and each signal processing operation that can be performed

Bandpass ac Filter and ae Deep I~-~1AmAmp Terminal Amplifier Isolation Amplif~ier PD'-8 C omputer- Controlled Com]ue n Hydrophones Marine Cable Attenuiator and C Peripea 1/3 octave Filter Equ imen I ~ ~~~~~~~~~~~~Bandpass Shallow Amp Terminal Amplifier Filteron L I~~~~~~~~solation Amplif ier (rJ Frequency I IFrequency Standard I 1 Synthesizer 1680 Hz Clock Signal Fig. 3.7. Equipment configuration at Bimini

Low Pass Samples of Filter Real Part | -- cosine phase Input l r 420 Hz | | Low Pass 420 Hz Osc illato Sampler, 210 samples/sec Cmplex Signal Sigr < sine phase __Low Pass | Samples of Filter - Imaginary Part Fig. 3.8. Complex demodulation of the hydrophone reception

39 at bandpass can also be performed on the lowpass complex signal with comparable results. (Van Trees [5] gives a proof of this duality for both deterministic and stochastic signals.) A distinct advantage of the complex representation in this case is that the number of samples per second which the computer must process depends only upon the bandwidth of the signal and not upon its center~ frequency.. After demodulation, the reception is passed through a number of digital filters, both to enhance signal-tonoise ratio and to perform the data analyses mentioned above; see Fig. 3.9. The transfer functions of the filters extract three distinct spectral components from the reception. The filter of Fig. 3.10a, called the carrier filter, extracts the carrier component in the reception. The bandwidth of this filter is sufficiently narrow to reject the surface modulation sidebands which surround the carrier, and yet is sufficiently wide to pass the long term fluctuations in carrier amplitude and phase which are of interest here. The filter of Fig. 3.10b is called the sideband filter; it has maximal response at frequencies corresponding to the transmitted lines in the probe signal spectrum. The

Complex Lowpass Signal Bandpass Carrier carrier power Reception XL -ilter and phase Carrier representative Removal 1 signal period Power Sideband sideband power -- Measurementd power Filter (Integrator) Power Noise Measurement noise power Filter (Integrator) to other filters Fig. 3.9. Block diagram of the signal flow for on-line processing

A k -.013 Hz H (f) (a)'cf 420 Hz HS~~~~~f*f4-01 Hz E- Af.t3z = b Hz 1~ f J V~~ 5 1 H (f)~~~~~.013 Hz A, xf=6Hz=T f Af \ f + Af c (b) c ~.013 Hz A - f H z~= N(f) Af Af f f + c 2 c 2 (c) Fig. 3.10. Transfer functions of the digital filters used in on-line processing (a) Carrier filter (b) Sideband filter (c) Noise filter

42 carrier line is removed from the reception prior to input to this filter. The total energy passed by this filter is a measure of the power received from the sidebands of the probe signal over a bandwidth set by the one-third octave filter of Fig. 3.7; for this experiment, the bandwidth was approximately 105 Hz. Each of the identical regions of nonzero response in the sideband filter is called a tooth,- in analogy to the teeth of a comb. The bandwidth of each tooth obeys the same constraints as the bandwidth of the carrier filter. The filter of Figure 3.10c is the noise filter. Its teeth are interstitial with those of the sideband filter so that it-has maximal response at frequencies where there was no power in the transmitted probe signal. The energy passed by this filter is a measure of the ocean noise over a bandwidth set by the input anti-aliasing filter, again approximately 105 Hz. The filters of Fig. 3.10 are actually implemented in the time domain instead of the frequency domain. After each 76.8 seconds of reception, the filter outputs are sampled, yielding: * Carrier filter: real and imaginary components of the carrier reception, coherently averaged over the 76.8-second accumulation period. * Sideband filter: 252 complex-valued samples of a representative period of the reception, minus

43 carrier, over the 76.8-second accumulation period. Noise filter: the total energy in the output of the noise filter, obtained by summing the modulus squared of the noise waveform over the 76.8second accumulation period. After sampling, all filters are reinitialized for the next sample period. Thus successive filter samples are uncorrelated. An on-line display showed the carrier amplitude and phase, the total sideband energy (obtained by summing the modulus squared of the output of the sideband filter), and the total noise energy. Both the displayed values and the carrier filter and sideband filter outputs were saved on magnetic tape. With these outputs, a representative period r(t) of the received signal at time t can be reconstructed for subsequent processing. 3.5 Selection of Data to be Analyzed The receptions from the Bimini receiver are of greater interest than those of the 12-km receiver since they have passed through a longer distance underwater. The recorded data from the Bimini receiving site span a period of 19 days. The shallow hydrophone data from the period 1130 hours to 1900 hours local time, 25 November, 1970, which will be called the test period, were selected for analysis for the following reasons:

44 (1) Digital tape recordings of sideband filter and carrier filter outputs are also available for the test period from the 12-km hydrophone -(Point 4 on Fig. 3.1). (2) The receiving and transmitting frequency standards were closely matched in frequency during the test period. In earlier parts of the experiment they were inadvertently about two millihertz apart, sufficient difference to compromise the phase structure of the transfer function. (3) The sideband power measurements over the test period show two distinct modes of behavior, also observed during single frequency experiments of previous years. One mode shows low, constant power which appears to be the normal behavior encountered in the Straits. (Figs. 3.3 and 3.4) The other mode shows higher, varying power. This mode usually appears during the autumn when a storm of long duration passes through the area from the north. Ostensibly, these two modes of behavior indicate a change in the underwater propagation path structure between Miami and Bimini. Some investigators have hypothesized that the anomalous behavior is due to a subsurface duct initiated

45 by the interaction of the storm with the surface layers in the Straits. The results of the on-line signal and noise power analyses during the test period are shown in Figure 3.11. The sideband power, labeled S, and the carrier power, labeled C, clearly show the two distinct modes of behavior mentioned above. Before 1440 hours, the S-power is low and relatively constant; after that time, it rises rapidly and is much more variable. In general, the carrier level tends to follow the sideband level (to be expected if the ocean is all-pass). The substantial fluctuations in the carrier level are due to cancellation caused by the phase interference phenomena explored in Chapter II. The signal-plus-noise to noise ratio, in dB, of the C- and S-power traces is the difference -between each trace and the trace labeled N, which is the wideband noise power measurement. During the test period, a number of photographs were made of the envelope of the multipath response of the channel. Tracings of these photographs are reproduced in Figure 3.12. It is difficult to tell the change in propagation mode in these pictures. 3.6 Transfer Function Extraction Given the above data, the transfer function of the Miami —Bimini channel was extracted. The transfer function

150 140.,...__. 130 1 120 1C.. 1 15 1 C 0g 110 100 _ _____ N 80... 11 12 13 14 15 16 17 18 19 Time (hours) Fig. 3.11. On-line power measurements for 25 November 1970

47 1. H H j uo = |\1030 hrs. Time 1.2 sec 0 ~ a? i, On 131.0 hrs. 1.2 sec rO 1750 hrs. 1.2 sec Fig. 3.1_2. Measured multipath responses of the Miarni-Eimini channel, 25 November 1970

48 of a channel, S(f), is defined in the frequency domain as S(f) = R(f) / P(f) (3.3) where R(f) is the noise-free channel response to input probe signal P(f). In this study a noisy version of R(f) can be found by Fourier transforming the measured time response r(t) of the channel. P(f) is the known spectrum of a CM signal. Applying Eq. 3.3 directly would cause the noise power across frequency in S(f) to be nonuniform, an undesirable effect. The transfer function definition used in this study, modified not to alter the noise level across frequency, is I! S(f) i = || R(f) | (3.4a) arg (S(f)) = arg (R(F)) - arg (P(f)) (3.4b) The available channel time response data consists of: 274 time series, each the output of the sideband filter for 76.8 seconds of underwater sound reception. Each time series s(t), t = 1,... 274, contains 252 complex points s(t, i) i = 1 252. 274 complex-valued points c(t), t = 1,... 274, each giving the average carrier amplitude

49 and phase over 76.8 seconds of underwater sound reception. The sampled channel time response r(t) is r(t) = c(t) + s(t), t = 1,..., 274 (3.5) formed pointwise by r(t, i) = c(t) + s(t, i) i = 1,..., 252 (3.6) A (time series, carrier) pair was produced every 99.6 seconds. Thus the data span 274 times 99.6 seconds, or 7.58 hours. 3.6.1 Fourier Transformation. Each 252-point Sfilter output time series s(t) was discrete-Fouriertransformed accoring to the following definition: 252 2Tv 252 252 J 252(i-1) (k-1) S(k,t) = X s(t,i) e, k = 1,..., 252 i=l (3.7) This operation was implemented using a special fast Fourier transform (FFT) algorithm. Using a 252-point transform, the eigenfrequencies of the transformed data lie exactly

50 on the transmitted probe signal line frequencies. Modifying the number of points in the transform to another value (say, 256) to expedite the FFT algorithm was considered unacceptable because such a match between transform eigenfrequencies and probe signal frequencies would not occur. As expected, the frequency spectrum S(k, t) had an extremely low carrier value. That the carrier value was not zero can be attributed to round-off error in the computation. 3.6.2 Reinsertion of the Carrier. The large amplitude of the carrier in relation to the sideband lines S(k, t) invited scaling of the carrier prior to its reinsertion into S(k, t). The scale factor chosen was such as to transform the transmitted spectrum of Fig. 3.6 into a pure sin(x)/x curve, i.e., the carrier was attenuated sufficiently to fit the smooth sin(x)/x curve of the other lines in the transmitted spectrum. For the CM signal used, the ratio of transmitted carrier power to the carrier power in the equivalent (sin(x)/x)2 power spectrum is p = 1 + L ( ) 1(3.8) For the sequence used, p = 62.03125. Since S(k, t) is a complex spectral value and not a power spectrum, each carrier value c(t) was divided by /p before

51 insertion into S(0, t): S(0, t) = c(t) (3.9) 3.6.3 Restriction in Frequency. The sideband power measurements of Fig. 3.11 indicate that the signal-plusnoise to noise ratio from the sideband filter varied from a minimum of about 10 dB to a maximum of about 25 dB. The centermost lines in the transmitted probe signal spectrum contain the most power. Subject to cancellation due to multipath phenomena, the signal-plus-noise to noise ratio of the centermost lines in the reception is then greater -than or equal to that of the sideband filter output. Given these considerations, thirty lines to each side of the carrier, as well as the carrier itself, were selected for further analysis and display. The power transmitted in these lines varied from 1.008 to 2.166 times the average power in a line in the main lobe of the transmitted probe signal. The frequency restricted transfer function is S' (f, t) for t = 1,..., 274:

52 S'(f, t) = S(k, t), k = 223,..., 252, f = k - 223 = 0,..., 30 (3.10) S'(f, t) = S(k, t), k = 1,..., 31,i, f = k + 30 = 31,..., 61 This transformation puts the carrier at f = 31. The prime on S' will be dropped for the remainder of this report. The duration of the time series s(t, i) is 1.2 seconds, yielding an eigenfrequency spacing of 5/6 Hz. Since 61 lines were selected by Eq. 3.10, the spectra S(f, t) span 50 Hz, i.e., 25 Hz to either side of the carrier. This bandwidth is 0.17 octave. At this point, the subtraction of the transmitted probe signal phase, Eq. 3.4b, was performed. 3.6.4 Phase Delay Compensation. Arg(S(f, t)) was modified to compensate for linearly increasing phase with increasing frequency, a phenomenon present in any wideband signal (e.g., Eq. 2.4). An interactive computer graphics program was written and used to display the phase of S(f, t) for each t. Provision was included in the program for subtracting a line of constant slope from the phase angles before display: arg(S(f, t)) = arg(S(f,t)) - (M-(f - 1) - B) displayed for t fixed (3.11)

53 The slope M and intercept B were continuously variable through potentiometers. Single values of M and B were chosen empirically with the assistance of this program to make phase as a function of frequency as flat as possible over the entire function S(f, t). The values chosen are M = 0.45507 revo lution/eigenfrequency B = 0.58593 revolution New phase angle values were then computed by performing the subtraction of Eq. 3.11 on S(f, t) for t = 1,..., 274. This phase delay compensation corresponds to shifting the time origin of the reception by 0.55550 seconds and to redefining the phase origin when computing arg(S(f, t)) 3.7 Summary of the Properties of S(f, t) The above sequence of transformations yields samples of the time-varying ocean transfer function S(f, t) in a selected frequency band, as defined by Eq. 3.4. The deviation between S(f, t) and the transfer function of Eq. 3.3 is in amplitude only. The deviation is graphed in Figure 3.13. S(f, t) contains 274-61 = 16287 complex numbers representing samples of the ocean transfer function over a 50 Hz bandwidth centered about 420 Hz and over a 7.58-hour period. The data may be viewed either as 61 complex-valued time

54 series each having 274 samples - S(f, *), or as 274 complex-valued spectra each composed of 61 frequency points - S(*, t). If cancellation due to multipath is ignored, the signal-plus-noise to noise ratio of each of the points is at least 10 dB. Plots of both the time series and the spectra are available in [7]. t 3 rl 2 > 0 11 21 31 l1 51 61 eigen frequencyv Fig. 3.13. Deviation in dB between Eq. 3.3 and Eq. 3.4 as a function of frequency

CHAPTER IV DEVELOPMENT OF INFORMATION DISPLAYS Efficient, human-engineered display procedures are a prerequisite to the visual analysis of S(f, t). The display procedures were specifically developed and optimized for this study since previously published procedures were unsuitable in one or more ways. The display procedures used are (1) parallel projection with hidden line elimination, and (2) contour mapping. These two techniques have seen widespread application ([12], [13], [14]), but infrequently with large numbers of regularly spaced data points as in S(f, t) By convention, S(f, t) is equated to the three-space equation z = f(x, y). Frequency increases along the x-axis, time increases along the y-axis, and z is the transformation of the transfer function data. The coordinate pair (x, y) at which a sample of z is available will be called a base point. The entire set of (x, y, z) coordinates will be called a data surface. 4.1 A Brief Introduction to the Display Techniques 4.1.1 Parallel Projection. Perhaps the best introduction to the parallel projection display procedures is by example. Figure 4.1 shows a parallel projection of a data surface containing 2000 base points. The angle of rotation of the viewpoint about the z-axis was chosen to be 45 55

56 ~x Fig. 4.1. Example parallel projection of a data surface with 2000 base points

57 degrees and the angle of elevation of the viewpoint above the x-y plane was fixed at 30 degrees. The reader may want to compare Fig. 4.1 with the rendering of the same surface on page 27 of [15]. The viewer must specify a number of parameters in addition to the coordinates of the data surface; chief among these are: The angle of rotation of the viewer about the z axis. The angle of elevation of the viewer above the x-y plane. - The aspect ratio of the rectangular grid of base points. ~ A parameter to control the vertical height of the picture. A vector of switches which control the production of baselines and cross hatching and the disposition of error comments. Appendix B explains the use and operation of the parallel projection display procedures in more detail. 4.1.2. Contour Mapping. In the equation of a data surface z = f(x, y) (4.1) one may choose to fix z at some value, say z0, and solve for y as a function of x: y = g(x, Z0) (4.2)

58 A contour line is defined as a set of connected points in the plane z = z0 which satisfies Eq. 4.2. Since g may be multiple valued, a single value of z0 may generate more than one contour line. The set of all contour lines generated by a value of z0 is called the contour set of g at z = z0. A family of contour sets can be generated by a family of z0 values. The projection of a family of contour sets onto the x-y plane is termed a contour map of the data surface f Figure 4.2 shows a contour map produced by the contour Fig. 4.2. Example Contour Map of Fig. 4.1

59 mapping display procedures from the surface of Fig. 4.1. Contour labeling has not been invoked in Fig. 4.2, although a primitive labeling facility is available. In addition to the data surface coordinates, the viewer must specify the following parameters: The family of z-values at which contours are to be mapped. *A vector of switches which control the drawing of border lines and the production of contour labels. * If labeling is invoked, the x-y coordinates of the locations of the labels and a value for the minimum distance between labels. Appendix C explains the use and operation of the contour mapping display procedures in more detail. 4.2 Design Criteria and Constraints Several constraints on the display procedures arise because the human visual perception mechanism is used to analyze the transfer function data. Interestingly, material published in this area by computer knowledgeable authors (e.g., Martin [11]) rarely mentions the perceptual aspects of the pictures, but concentrates on the system aspects surrounding the displaying of pictures with few base points (typically <500). Psychological research [26] has shown that a human observer asked to do feature extraction on a complex visual presentation (such as a parallel projection of a data

60 surface) first builds an internal perceptual model of the objects in the presentation. The observer then extracts features from the model rather than from his retinal image. Several different visual cues are used to build the model. As texture and motion cues are the most powerful, particular attention should be paid to these cues when constructing pictures for human observers. Humans also appear to have neurophysical mechanisms which are especially adapted to discern discontinuities in a visual presentation. Human observers are thus adept at discerning peaks and edges in complex visual presentations. To take advantage of these mechanisms, data values of a priori interest to this study (e.g., amplitude maxima and minima) should be transformed into peaks and edges in the parallel projection pictures to be observed. Contour maps of the transformed data can also be used to track the path of the edges in the f-t plane. Also considered were the constraints imposed by characteristics of the available plotting devices. The devices fall into two categories: the interactive graphics terminal using either refresh or direct view storage tube technology and the digital incremental plotter (including COM - com — puter output on microfilm). The distinction here is between low resolution, small picture size, instant turnaround devices and high resolution, large picture size, long turnaround devices.

61 Under these constraints, the main issues in picture production are: * The costs involved * Size of the picture as it affects texture cues The aspect ratio, orientation, and rendering of the picture on the plotting device surface, as they affect texture and motion cues * For the parallel projection, the variation of the picture in the vertical direction Use of an off-line, high-resolution digital plotter in conjunction with an on-line, interactive graphics terminal addresses all these issues. Each is discussed in detail below. One component of the cost question is the price of generating the pictures themselves. Initially, the investigator will not know the optimal presentation of the data. He must iteratively adjust various picture parameters to highlight features of interest and to produce an informative picture. If the cost of producing a picture is excessive, he may choose to terminate his investigation prematurely, possibly missing an insight into the data. Consequently, there is a strong incentive to produce fast display procedures. This goal is especially important due to the large number of base points in S(f, t). Previously published display procedures could not be easily adapted to this goal. For information on the software engineering

62 methods used to produce the display procedures, see Appendix A. For the investigator utilizing an interactive system, another component of the cost is the charge for connect time to the computer. There are periods during which the computer is idle while the investigator is contemplating a picture. Relatively high connect charges may also lead to the premature termination of the investigation. An off-line, digital plotter combined with fast display procedures answers the cost issue. The plotter is used to produce a hard copy of those pictures which the investigator deems sufficiently interesting when browsing in the on-line, interactive part of the investigation. While online, he interactively adjusts various picture parameters, thus assuring quality hard-copy pictures for more intensive, off-line study. Concern with picture size is justified by the large number of base points in each picture. Plotting a 16,000base point picture on a plotting device which only has 1024 addressable points along each axis can lead to unacceptable distortion. Moreover, the lines drawn by any plotting device are not infinitely thin. If, as in typical interactive graphic terminals, the horizontal axis is less than, say, ten inches long, texture cues may be lost due to blurring and overwriting of the lines comprising the picture. To overcome these technological limitations, the

63 viewer would ideally like both to inspect small sections of the picture, magnified to full-screen size, and to view the whole 16,000-point picture as an entity. Although the technology to expand sections of a picture (viewporting) is well-understood and commonly available on interactive terminals [12], large screen interactive terminals capable of faithfully portraying complete 16,000-point pictures are rather expensive. Again, the digital plotter is of use. The parallel projections drawn by the plottermeasure approximately 24 by 29 inches. The contour maps from the plotter measure approximately 12 by 28 inches. These sizes are large enough to give high resolution (5000 addressable points along each axis) and to counteract the thickness of the plotter pen (0.5 mm). The pictures are amenable both to viewing of local regions and to viewing of entire pictures as single entities. Special reproduction methods have been used for the reduced size pictures in this report so as not to lose the viewing properties of the originals, The viewer must be able to alter the aspect ratio of the picture at will. A picture which is too wide or too narrow can restrict the viewer's perception either by occupying too large a visual field or through loss of texture. In a parallel projection, a desirable feature is the capability to rotate the data surface in three dimensions.

64 The viewer may thus explore the shape of various features in the surface by viewing them from different directions. Ideally, this rotation should be done dynamically in real time [26]. Unfortunately, there are so many base points in S(f, t) that dynamic rotation is effectively impossible. The burden of evoking three-dimensionality must perforce fall to texture cues in static pictures of rotated data surfaces. Some choices of rotation angles and aspect ratios can lead to poor pictures. Figure 4.3 shows five views at various rotation and elevation angles of a data surface developed for this study. The rotation and elevation angles- are indicated alongside each view. At low angles, the lines comprising the picture lie quite close to each other; consequently the width of the line obscures features. (This effect could be compensated somewhat by choosing the aspect ratio to make the picture longer.) In addition, at low elevation angles, features close to the viewer obscure features farther away. At high angles, many features in the data are lost. Air travelers notice a similar effect when flying over a mountain range —it looks far less rough than it does from the ground. The best elevation angles seem to lie between 20 and 40 degrees. Within this range the exact choice of angle makes little difference. The choice of rotation angle depends upon the rendering of the data surface. When the surface is rendered by lines parallel to the x axis,

65 (20,30) (10,15) (45,60) Fig. 4.3, Five views of the same surface at different rotation and elevation angles. The angles are marked next to each view in the form (rotation, elevation). Because of the halftone reproduction process used, these-views are more pleasing to the eye than the originals

66 the rotation angle is best kept below 45 degrees. When the surface is rendered by lines parallel to the y axis, the rotation angle is best kept above 45 degrees. In the vicinity of 45 degrees and for pictures which are not visually crowded, rendering by both x and y lines (cross hatching) is most informative (c.f., Fig. 4.1). Figure 4.4 illustrates the deleterious effect of improper choice of rendering. The data surface is rough and peaky (Fig. 4.4a), and the location of the peaks along the y direction is of interest. Cross-hatching this surface (Fig. 4.4c) darkens the picture and obscures many of the smaller peaks by destroying texture cues. This rendering effect is not so pronounced in smoother data surfaces, but is still present (See Fig. 4.5.). In Fig. 4.4b, the lines in the picture parallel to the y axis are quite superfluous to the rendering. They were included simply to indicate equal divisions along the x axis. Yet these lines capture the eye of the observer and make it difficult to follow the peaks in the y direction. Plainly, aspect ratios, rotation and elevation angles, and rendering methods should be carefully chosen to elicit informative pictures. In a parallel projection, let the "vertical variation" be defined as the maximum z value in the data surface presented for display minus the minimum z value in the data surface presented for display. The value of this

(a) (b) (c) Fig. 4.4. Three views of the same data surface showing the effects of different renderings. (a) No crosshatching (b) With frequency markers (c) With crosshatching

Fig. 4.5. Two views of a relatively smooth data surface, showing the effects of different renderings

69 variation measures the vertical amplitude in the picture. The vertical variation in a picture significantly influences perception of features. If the vertical variation is quite small, the surface appears too smooth, and extracting features is correspondingly made difficult. This is the problem in the (45, 60) view in Fig. 4.3. Conversely, if the vertical variation is too high, the picture will appear noisy, obscuring significant features. Figure 4.6 illustrates this point with three pictures of the same data surface but with increasing vertical variation. The vertical variation of a picture must be adjustable to allow the most informative presentation. 4.3 Conventions for the Parallel Projection Pictures Since certain choices of rotation and elevation angles and grid aspect ratio were found to produce consistently pleasing pictures and since it was desired to have some uniformity of presentation to facilitate inter-picture comparison, nearly all the parallel projections in Chapters V, VI, and VII are oriented as shown in Fig. 4.7. The angles of rotation of the pictures in Chapters V and VI is usually 20 degrees, while the rotation angle for the pictures of Chapter VII is 15 degrees. For all pictures, the angle of elevation is 30 degrees. The aspect ratio of the base point grid is one to one. No single choice of the vertical variation was suitable for all the pictures. The vertical variation in the

Fig. 4.6. Three views of the same data surface, with increasing vertical variation from left to right

71 18::53 Hrs 395 Hz yuen 11:18 Hrs 445 Hz Fig. 4.7. Example transfer function parallel projection, with coordinate system conventions

72 pictures ranges from 8.75 per cent to 25 per cent of the length of the x axis. Others may find a different choice of parameters more pleasing. The point is that the option to change the picture parameters is readily exploited if desired. The coordinate system for the parallel projections is also shown in Fig. 4.7. The original x axis is frequency in Hertz, the original y axis is time, and the z axis is whatever function is captioned under the picture. The 420 Hz carrier is centered on the x axis, 4.4 Conventions for the Contour Maps All the contour maps in Chapters V and VI are oriented as shown in Fig. 4.8. The x axis is frequency in Hertz; each tick on the x axis denotes five eigenfrequencies — 4+6 Hz. The carrier frequency is centered on the x axis. The y axis is time. Each tick on this axis corresponds to 8.3 minutes (exactly). The outermost contours surrounding a region are labeled with the value of the z coordinate of the data surface named in the picture caption. To select z coordinate values for the contours, a program was written to scan a data surface and print its maximum value, its minimum values, and its dynamic range in dB. Contouring values were selected on the basis of this information (tempered by physical considerations) to evenly straddle the z coordinate range in the data surface. In a few cases, z values toward one end of the range were

73 18:-53 Hrs LI ~ g la s3 *L II: o. W':,'5' -h. T-rX- 11: 18 Hrs 395 Hz frequency w 445 Hz with coordinate system conventions with coordinate system conventions

74 selected to emphasize features of interest.

CHAPTER V ANALYSIS OF THE TRANSFER FUNCTION AMPLITUDE Experimental observations of both the amplitude and phase of a single-frequency reception have been made for many years in the Miami-Bimini range [10]. The magnitude and timing of the amplitude variations at a single frequency do not exhibit an obviously consistent underlying behavior. However, as shown below, when viewed over a wide passband instead of a single frequency, the evidence is strong that the amplitude minima and maxima behave as predicted by Eq. 2.15. The method used is to apply a tranformation to the amplitude II S(f, t) to emphasize regions of interest by turning them into peaks and edges. The regions are then studied visually by viewing parallel projections and contour maps of the transformed data. 5.1 Amplitude Maxima The square of the amplitude was chosen to study amplitude maxima. This transformation was chosen because it emphasizes higher amplitudes at the expense of the lower and because it computes a statistic of intrinsic interest: the time-varying power spectrum. The time-varying power spectrum, | S(f, t) | 2 is shown in Figure 5.1. Features along the time direction tend to follow the sideband power observation of Fig. 3.11. 75

~Tj tq Ul H-3 (D 0 cD (D C Q'Cl rt tI —, r-j v ~ ~ ~ ~ ~ ~ ~ ~ ~ VI/r/l///l////~)))) 9 r(l(llllll(llllll, ((///1/1'//I((~Kl('u~

77 The time of the propagation mode change mentioned in Chapter III is shown by the arrow at the right of Fig. 5.1. At this time there is a small upward step in amplitude which might be overlooked at first glance. For approximately 45 minutes after the event, the power spectrum qualitatively remains unchanged. Plainly, the power spectrum is not a key indicator of propagation mode changes. For the first half of the test period (including the 45 minutes after the propagation mode change), the power spectrum is low and relatively smooth. Some frequencies are received somewhat more strongly than others, but the variation across frequency is not striking. At the beginning of the second half of the test period, the amplitude rises sharply and becomes strongly differentiated across frequency. For a time there are three strong frequency bands, five to ten Hz wide. These rearrange into two strong bands and then die out. The total time duration of these strong receptions is a little over two hours. The remaining portion of the power spectrum shows behavior akin to the first half except that there is a noticeably strong but diffuse band of receptions centered about twelve Hz below the carrier. One feature of interest is the consistent amplitude response at 436 Hz, indicated by the arrow at the bottom of the figure. There is significant amplitude at this frequency nearly 80% of the time. Also of note is the consistent lack of power at 439i Hz, just four

78 eligenfrequrencies, above 436 Hz. Even. during time periods of high amplitude response, there is very little power atthis frequency. The highest amplitude reception occurred on the. carrier line about halfway through the test period. Let this. amplitude be called SMAX. From Fig.; 3.11, the sideband filt.er output signal+noise/noise ratio at. this point was approximately 28 dB. This highest amplitude point is the "moun-tain: peak" near the center of Fig. 3.1, at f = 31 and t = 1.57 Figu~re 5.2 shows contours of the power spectrum in units of decibels down from SMAX. Contours are located at 5, 10, and: 15. dB down;: all the-se levels are at least 10 dB above the noise level in Fig. 3.11.. Patterns are not as easily seen in the contour map as in the parallel projection. However, close inspection shows sporadic movement of amplitude contours along diagonals in the f-t plane. Such movement is confirming evidence for Eq. 2.15; spectral amplitude at time t is reproduced at t + At bult shifted in frequency. 5.2 Amplitude Minima The track of amplitude nulls in the f-t plane offers even stronger confirmation for Eq. 2.15 since amplitude nulls are much narrower in bandwidth than in amplitude maxima. One possible transformation to turn amplitude nulls into peaks is relative path loss statistic L:

79 0 1~ sW/t is. I. FE <,., A., E Fig. 5. 2. Contours of the time-varying power spectrum at 5, 10, and 15 dB down from SlAX 4 Q) Y

80 L(f, t) = -20 logo10 S(f, t)J /SMAX) (5.1) L(f, t) is a measure in dB of the relative loss introduced by the ocean at frequency f and time t. It emphasizes the lower amplitude portions of the transfer functions at the expense of the higher. A parallel projection of L(f, t) is presented in Figure 5.3. The. movement of nulls along diagonals in the f-t plane is quite apparent. The nulls sweep toward lower frequency in the front portion of the data surface. In the middle portion of the data, there is a "horseshoeshaped" sweep of the nulls, where the nulls converge toward the center of the band with increasing time. At the rear of the surface, there appears to be some motion of the nulls toward the right. Plainly, inspection of adjacent frequencies allows an understanding of amplitude nulls in terms of their diagonal movements in the f-t plane. Such an understanding is impossible to glean from singlefrequency data. Because they are clustered in the f-t plane, a contour map can give a concise picture of the track of the nulls. Figure 5.4 shows contours of L(f, t) at 20, 25, 30, and 35 dB loss. There is no obvious pattern in the location of the nulls per se. However, the movement of the nulls along diagonals in the f-t plane is quite apparent.

81S Fig. 5.3. The relative path loss function

-,-I 0 4Jr< 40-~~~~~~~~~~~c cc) 00 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-4 tIII t. I,. rL-~~~~~~~~ 0' ~~~~~~~~~ >rd~~~~~ C~* ~" Co~ ~ ~ ~~~~~~~%**o*-'..)c o'o *~~~~~~~~~~~~~~~~ oLC) -H1

83 5.3 Very Low Amplitudes The average level of the ambient noise, estimated from Fig. 3.11, lies on the 25 dB contour in Fig. 5.4. Eliminating the 20 dB contour yields Fig. 5.5 whose outside contours are at 0 dB signal-to-noise ratio. Inside the enclosed regions in Fig. 5.5, the nulls go deeper than the 25 dB ambient noise level. The highest value attained by the relative path loss is 59.6 dB; thus there is a value in S(f, t) which is 34 dB below the estimated wideband noise level. The ambient noise level should serve as a ceiling on the value of L(f, t) since the S-filter of Section 3.4 responded to noise as well as signal. Below a signal-tonoise ratio of 0 dB, points in S(f, t) are essentially samples of the ambient noise. Assuming the noise is a Gaussian random process, II S(f, t) jj 2 will have the Chisquared distribution with two degrees of freedom below 0 dB S/N. Letting P = I S(f, t) 2, the probability density of P is fp(x x 1 2 fp(x) = e, 0 x < (5.2) so that p p _P_ Prob(P<p) = fp(x)dx = 1 - e 2 (5.3) 0 To

TGAGa aSTOU pupq -apTM pGqwuITs tai O14 4 ST Ino:uoo gp SZ aL, *sso{;qpd GATqtTEl gp SE pup'OE'JSZ q uoIToung ssoT qgcId OATeiaa aq Jgo sO noS uoOD *S*S *bT 0-sa ~0s oo 0 0;: sz X 0 <) t )0-.s0 d' _ _ ~ _' " t4 ob,

85 The expected value of the noise power is 00 p P 2 E(P) = ] e dP = 2 (5.4) 0 and Prob(P<2) = 1 - - = 0.632... (5.5) e There are 1303 values in S(f, t) with amplitude at or below the noise level. If we assume that these 1303 values comprise the 63.2% (Eq. 5.5) of the population at or below the average noise power, then the total population size is 2061 samples. The solid line in Fig. 5.6 is the expected variation of the number of values in the population below a given noise power level. It is the product of the population size and the Prob(P<p) from Eq. 5.3. The measured values given in Fig. 5.6 are quite close to the expected values. This is substantial evidence to indicate that at very low amplitudes, S(f, t) is driven by Gaussian noise. 5.4 Conclusions Parallel projections and contour maps of transformations of the wide bandpass transfer function amplitude have yielded the following insights into the time and frequency variations of the transfer function amplitude.

86 2000 0 1000 _ 500 100 U) _ 0 0 - 4-4 1100 o - = expected k = _ O = measured 50 rd 0)/ ~ d 20 0 2 0 / -40 -30 — 20 -10 -0 Amplitude in dB relative to the wideband noise level Fig. 5.6. Expected and measured cumulative distribution functions for those amplitudes below the estimated wideband noise level

87 ~ Minima and maxima in the transfer function amplitude often move on diagonals in the f-t plane, as predicted by Eq. 2.15. * The behavior of S(f, t) cannot be fully explained by Eq. 2.15 as there exist frequencies at which the path loss is consistently low (e.g., 436 Hz) and frequencies at which the path loss is consistently high (e.g., 439 Hz). Changes in the mode of underwater sound propagation are difficult to observe solely on the basis of variations in the amplitude of the transfer function. The ambient noise drives the transfer function amplitude when the signal-to-noise ratio drops below 0 dB. Moreover, a Gaussian model for the noise is strongly supported by the behavior of the amplitude below 0 dB signal-to-noise ratio. The following conclusions concern the graphical analysis techniques themselves: Parallel projections of the time-varying power spectrum are useful in studying the high-amplitude features in the transfer function. Contour maps are of marginal utility in studying the high amplitude regions of the time-varying power spectrum. The path loss transformation adequately emphasizes interesting low-amplitude features in the transfer

88 function. Both contour maps and parallel projections are useful techniques to study the path loss, and thereby to study low-amplitude regions of the transfer function.

CHAPTER VI ANTALYSIS OF TIHE TRANSFER FUNCTION PHASE Observations at a single frequency have shown that the phase of underwater sound receptions usually varies quite slowly in time. The exception occurs when the re — ceived amplitude drops so far into the ambient noise that reliable phase measurements cannot be made. In this case, the phase measurements before and after the amplitude null often differ by 90 or 180 degrees. These jump discontin — uities are usually imputed to changes in the path struc — ture between source and receiver, but since noise levels were not measured in these observations, some uncertainty remains as to the cause. As shown below, when viewed over a wide passband instead of a single frequency, the phase exhibits a high degree of structure in both the frequency and the time directions. Anomalies due to low signal-to-noise ratio become quite distinguishable from variations due to path structure changes. In this study, phase angles will be measured in revolutions (rev.) as opposed to radians or degrees. For example, angles ranging from 0 to 540 degrees will be said to range from 0 to 1.5 revolutions. The visual analysis methods of Chapter'V must be extended to cope with the phase, as amplitude measurements 89

90 are continuous by nature, but phase measurements are not. The phase of S(f, t) is computed via -1 Im(S(f, t)) f(f, t) = tan t)) (6.1) Since the arctangent function is multiple valued, a convention must hbe made regarding principal values in the com — putation of Eq. 6.1. Figure 6.1 shows ~(f, t) where principal values lie between 0 and 1 rev., and Fig. 6.2 shows f(ff t) ranging from -.5 to +.5 rev. The raw data comprising the pictures are the same, yet the perceived features are quite different since the principal value convention radically changes the location of t'he peaks and edges in the pictures. The problem is then to define a transformation on f(f, t) which will produce informative pictures and which is invariant under the principal value convention. 6.1 Reconstructed Phase Angle From purely physical considerations, the phase of S(f, t) must vary continuously over a multi-revolution range as the sound ray paths from source to receiver change in length (Eq. 2.4). Thus the computed phase of S(f, t) should extend over the full Rierann Surface of the arc — tangent function. Ficure 6.3 shows this surface with the principal value of the arctangent chosen from —.5 to +.5 rev.

91 Fig. 6.1. The raw phase p(f, t) with principal value convention 0 to 1 revolution

92 Fig. 6.2. The raw phase ~(f, t) with principal value convention -.5 to +.5 revolution

93 chase angle the oath of - 1 revolution Fig. 6.3. The Riemann Surface of the arctangent function

94 Under strong signal-to-noise ratio conditions, the phase at a fixed frequency is observed to change at most about 0.1 rev./minute [10]. Thus it should be possible to track the phase, reconstructing its movement over the full Riemann Surface from frequency samples in a 1 rev. range. The reconstruction algorithm used to track the phase is explained in Appendix D. It reconstructs the phase at a selected frequency independently of the behavior of S(f, t) at any other frequency. The algorithm was applied to t(f, t) to yield the reconstructed phase 6(f, t) whose value ranges over the real line. However; the initial results were mixed (see Fig. 6.4). Since the initial conditions for the algorithm were chosen arbitrarily to be zero, there were notable variations in phase across frequency. Such variations were not present in the pictures of spectra used in Section 3.6.4 to determine M and B for Eq. 3.11. To correct for the across-frequency phase variations in Fig. 6.4, S(f, t) before and after the propagation mode change was treated as two separate entities. The two spec — tra S(f, 50) and S(f, 158) were determined to have acceptably high signal-to-noise ratio at all frequencies. Correction functions C(f, 50) and C(f, 158) which were constrained to operate only on the integer part of 0 were empirically determined to minimize the across.-frequency variation. The correction functions were then applied to all of 8'

95 Fig. 6.4. The phase e(f, t) as reconstructed by Birdsall's algorithm

96 0(f, t) = e(f, t) + C(f, 50), t = 1,...,103 (6.2) 0(f, t) = 0(f, t) + C(f, 158), t = 109,..,274 (6.3) These corrections could also have been made by choosing proper initial conditions for the reconstruction algorithm applied independently to the two segments of f(f, t) Figure 6.5 is a picture of the corrected 0. The time of the propagation mode change is marked by the arrow at the right of the figure. The change is readily apparent as a feature in the picture due to the jump discontinuity in ~ at t = 109. The mode change is much more apparent in the phase picture than in any of the amplitude pictures. Contours of Figure 6.5 are presented in Fig. 6.6. The arrow at the right marks the propagation change. Due to the jump discontinuities in j r contours before and after the change are not colinear. The striking feature in Fig. 6.6 is that the contours move coherently across frequency. Immediately before the propagation mode change, they are moving strongly toward higher frequencies, a movement continued after the change. Toward the end of the picture, the upward movement ceases and becomes a wseak downward movement. This coherent movement of the contours in the f-t plane is supporting evi — dence for Eq. 2.15.

97 g Fig. 6.5- The corrected phase 0 (f, t)

(;'j)9 O ssLqd paoDosloD e14; _o sjnouo3'*9'9'frc *I I O ~~~~~~~~~~~~~~~~~~~~d-, 0. - 0 o' o~ o o 9- 0 tI' 1f[0 *' 0. _~~._ 0 0 O o ( o rJ 1_ ~ 86 -~~~~~~~~~~~~~~ i _ _. ~ 111 -tt o 0"~~~~~~~~~ P II~ ~ ~ ~ ~ ~ o C-. o',6' - t s; O II ii" "~~~~~~~~~~~~~~~~~~<u _~~J P O86

99 The sequence of transformations from Fig. 6.1 to Fig. 6.5 preserved the fractional portion of each phase measurement. The only liberty taken was to move a measurement from one sheet of the Riemann Surface to another. It is remarkable that the transfer function phase is so well structured that such simple transformations will cause the phase surface to appear so flat and the contours to behave so coherently. 6.2 Problems in Corrected Phase Due to Noise A significant number of the features present in Figs. 6.5 and 6.6 are due to low signal-to-noise ratio in S(f, t). As an example, consider the region bounded by t = 155 to.170 and f = 36 to 38. Tables of the path loss and corrected angle for this region are shown in Fig. 6.7; a plot of the time track of S(f, t) in the complex plane is shown in Fig. 6.8. Frequency line f = 37, the line of interest, is indicated on both Fig. 6.5 and 6.6 by an arrow at the rear of the picture. In Fig. 6.5, the selected time interval corresponds to the region where the phase rises from the plane into the ridge which persists through the remainder of the picture. Observe the corrected angle table for f = 37. The values traverse a one-revolution range. The adjacent frequencies do not traverse this range either before or after the traversal in line 37. This is a seeming violation of Eq. 2.15.

Path Loss (dB) & (rev.) f + f + 36 37 38 36 37 38 155 7.7 14.0 15.2 155 -.8.2 -.9 156 5.7 13.1 11.8 156 -.8.1 -.9 157 5.6 12.2 9.4 157 -.8.2 -.9 158 7.3 19.6 9.5 158 -.9.1 -1.0 159 7.5 21.8 8.1 159 -.9.0 -1.0 160 7.7 23.7 8.1 160 -.9.0 -1.0 161 9.8 36.0 8.9 r 161 -.9.3 -1.0 ~ 162 9.7 40.9 9.6 ~ 162 -.9.6 -1.0 163 9.9 32.6 10.1 163 -.9.8 -1.0 164 9.6 24.8 9.0 164 -.9.9 -1.0 165 9.9 29.7 10.6 165 -.9 1.1 -1.0 166 10.2 22.7 9.7 166 -.9 1.0 -1.0 167 10.1 15.5 9.8 167 -.9 1.0 -.9 168 12.3 15.3 9.9 168 —.9.9 -.9 169 13.2 12.3 8.5 169 -.9.8 -1.0 170 10.2 10.9 7.5 170 -.9.9 -.9 Fig. 6.7. Path loss and phase angle of S(f, t) during a spot fade on line f = 37

OLI oq SSI =4 Jog auPid xa{lduIoD aOq uT 8 puP'LE'9~ = 4.o4 (J')gS _o o 4' -8*9 *BTJ L L 9E NI, TOT L -,-4-,,,- 9~~~~~~~~~l

102 The key to the behavior of line 37 lies in the path loss. Recall that the ambient noise lies at 25 dB path loss. Figure 6.7 shows the signal-to-noise ratio on line 37 dropped well below 0 dB, and thus the measured phase angles are essentially the phase angles of noise alone. That they ascend in so monotonic a pattern must be ascribed to chance. Yet this ascent introduces an eye-catching feature into Fig. 6.5 by introducing a ridge lasting over the rear half of the picture. In Appendix E, the probability density function f(q) for the phase Q of a signal in additive Gaussian noise is derived, conditioned on the voltage signal-to-noise ratio v. Graphs of f(P) are shown in Fig. 6.9 for various signal-to-noise ratios from -10 dB to +10 dB. As the signal-to-noise ratio decreases, the measured phase has increasingly wider variance. A graph of 90% confidence intervals for measured phase versus signal-tonoise ratio is given in Fig. 6.10. The 90% confidence interval is about 80 degrees wide when v = +5 dB. Consequently, the outer contours in Fig. 5.4 can be interpreted as +40-degree confidence contours for the measurement of phase. In like manner, the outside contours in Fig. 5.5 are +80-degree confidence contours on phase. 6.3 Editing 0 for Minimal Variation Using the contour maps of path loss in Chapter V as phase confidence maps, the regions of the f-t plane with

+10 dB l.75 1.5 0 +6 dB -P +3 dB r-+ 0 6dB 0.75 rb I nl \I IL 0 -10 dB 0.25 -160 -120 -80 -40 0 40 80 120 160 Phase angle in degrees Fig. 6.9. The probability density function of phase angle as a function of signal-to-noise ratio in dB

+10 ro + 5 0 0 10 0~~~ l I I!! I I I I I I I I I I I I! I I 4-J -5 -160 -120 -80 -40 0 40 80 120 160 Phase angle in degrees + Fig. 6.10. 90% confidence limits on measured phase as a function of signal-to-noise ratio in dB

105 phase confidence intervals of 80 degrees or larger were marked on a table of the e(f,. t) function. The integer part of each 0 measurement in these regions was altered manually to minimize the variation in both the t and the f directions. No algorithm or measure of smoothness was developed to evaluate the alterations. Instead, 0 was viewed in parallel projection as editing proceeded to evaluate the effect of alterations. Figure 6.11 shows the & function after editing. Several features due to noise which are present in Fig. 6.5 are missing from Fig. 6.11. The contour map of the edited 6, shown in Fig. 6.12, is much smoother than the map of Fig. 6.6. The contours are not interrupted so often by noise-induced features along the t direction. Although human editing has improved the 6 pictures, it cannot be recommended as it is a tedious process subject to error and is not objective. A preferable procedure is to define a set function which can be applied to the whole of 6 without iterative selection of the desired sheet of the Riemann Surface in the low signal-to-noise ratio regions. The selected function should also, of course, transform informative regions of 6 into peaks and edges. 6.4 4The Phase Rate There are several desiderata for a transformation of 0 which is informative and yet insensitive to an incorrect choice of sheet in the Riemann Surface of phase.

106 Fig. 6.11. The phase 8 (f, t) after manual editing to remove features due to noise

(;'3)9 esnqd pe4TpG eq 3o slnouoD T g~ ~ O 00 -Io s ~SK ~ ~ ~s't. s't'4. - Os- s Is' 7 ~~~~~~~~~~~~~~~~~'l Q~~~ BQ 2 0 I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0 0 4f S~ So~ -S S~b ~ _ 0~ _ _,o _.1, /~~~~~~~~. O _ o a~~r O~ ~' o o~lOT'

108 Assignment of 0(f, t) to the wrong sheet introduces a step comrponent into ()(f, t). Consequently, the selected transformation should be insensitive to the mean level in G(f, t) * Since the object of this study is in part to exhibit the stabilities and instabilities in underwater sound propagation, the selected transformation should clearly distinguish between areas in the f-t plane where rapid changes are taking place and areas where gradual changes are taking place. * Past measurements of the frequency content of the phase of single-frequency receptions Show that the power spectrum of the phase falls rapidly with increasing frequency, in a manner reminiscent of the rectangular hyperbola. To compensate for this decline, the selected transformation should emphasize the higher frequency content of 0 at the expense of the lower. The selected transformation computes the phase rate, 6(f, t), so named because is is the slope of the linear regression line between several consecutive values of 0 and time, at fixed frequency. More formally, n consecutive values of 0(f, t) were regressed against time t to yield:

109 [nk-6 (f,k+t-l)]- [-k] [[0(f,k+t-l)] (6.4) [nXk2] - [k] 2 t = 1,..., 275-n where all summations run from k = 1 to n Equation 6.4 fulfills the three goals above. Plainly 6 will have higher amplitude where 6 is changing rapidly than where it is changing slowly. To see that the other two goals are satisfied, observe that Eq. 6.4 conducts a filtering operation on 6 to produce 0. An expression for the magnitude of the transfer function || H(f)|[ of the phase rate filter of Eq. 6.4 is derived in Appendix F and graphed in Fig. 6.13 for n = 8. There is a hard zero in H(f) at zero frequency and thus the phase rate is completely insensitive to the mean value of 6. H(f) also climbs at 6 dB per octave from zero frequency and passes high frequencies more readily than low frequencies. The value n = 8 for Eq. 6.4 was chosen because the major features in Fig. 6.11 require about eight time steps to change character. This choice was one of convenience rather than necessity and should be viewed as subject to change if warranted. The phase rate filter was applied to the corrected 6 of Section 6.2 to yield Fig. 6.14. These pictures show the phase rate before the 6 editing of Section 6.3. Consequently, some of the features are due to noise, but the

110 0 - (a) 25 20 folding 50l~~~ -*~~~~~~~~~~~~~~~ _ / \frequency,0 cn 30 20 15 12 10 9 8 7 6 5 4 3. 32 120 Period in Minutes n- (b) folding 0I frequenc~.r) 600 300 2V0 180 120 90 G0 30 20 15 12 1 0 98 7 6 5 4 3.32.Duration of Period in Mlinutes plot () Frequency-period conversion nomogram " n co a! 60 C C)'o C- C C'~'" -4,, Frequency in MilliHertz

111 (a) Fig. 6.14. The phase rate O(f, t) of the corrected 8 (before editing). (a) At 20~ rotation (b) At 0~ rotation (Continued on next page)

112 (b) Fig. 6.14. (Continued)

113 effect of the noise is localized. For example, observe the region where the deep fade occurred on line 37. e rises noticably here, but falls rapidly after the fade has passed. Features in e which persist over long times become localized in 0 Applying the phase rate filter to the edited 0 of Section 6.3 yields Fig. 6.15. In these pictures, features in 6 due to noise have been largely removed. The propagation mode change is immediately apparent in the pictures, extending across the whole frequency band. e is moving toward lower frequencies before the propagation mode change (more easily seen in Fig. 6.15b), and is moving toward higher frequencies at the end of the test period. The horseshoe feature of Fig. 5.3 is also easily seen. An unexpected feature is the uniform undulation in 6 parallel to the frequency axis throughout the whole picture. A measurement from Fig. 6.15 shows the period of these waves to be about 16 minutes. Consequently, they cannot be associated with the choice of n as 8 in Eq. 6.4 and must be present in the data. It is significant that these 16minute waves do not appear in pictures of 6 (e.g., Fig. 6.11), but are so clearly present when the high-frequency emphasis of Eq. 6.4 is applied. One possible oceanographic explanation for the 16minute waves in e is that they are an effect of internal waves in the Straits. Since the phase is a measure of the travel time between source and receiver, the phase rate

114 (a) Fig. 6.15. The phase rate G(f, t) of the edited 8 (a) At 20 rotation (b) At 0~ rotation (Continued on next page)

115 (b) Fig. 6.15. (Continued)

116 measures variations in that travel time. Internal waves, being small, periodic perturbations in the speed of sound, cause small, frequency-independent variations in travel time. These in turn should appear as periodic, frequencyindependent variations in the phase rate. The 16-minute waves are studied in rmore detail in Section 6.6. To handle the remaining features due to noise in Fig. 6.15, observe that only a small portion of the f-t plane lies within the low signal-to-noise ratio contours in Fig. 5.5. Since the effects of noise are localized in 0 few features can be changed if 0 is simply set to zero in these regions. In Fig. 6.16 e from Fig. 6.14 was set to zero whenever the estimated signal-to —noise ratio fell below +3 dB. Comparing these two figures shows immediately which features are associated with low signal-to-noise ratio (diagonal movements of 6) and which are independent of signal-to-noise ratio (the propagation mode change jump and the 16-minute waves). 6.5 The Power Spectrum of the Phase Rate The two-dimensional power spectrum of 0 will show the presence in the 6 pictures of ridges which follow diagonals in the f —t plane. The power spectrum was determined by computing the two-dimensional Fourier transform of 6. To speed the transformation process, a subset of the 6 picture in the region t = 1 to 256 and f = 1 to 60 was transformed.

117 Fig. 6.16. The phase rate obtained by setting e(f, t) from Fig. 6.14 to zero wherever the signalto-noise ratio fell below +3 dB

118 The power spectrum of the 0 was computed using both coherent and incoherent averaging. Define the function F as follows: 60 tmax F(u, v, tofs, t) = (f, t+ tf) ofs' max ofs f=l t=l (f-l)(u-l) + (t+tofs e71 G 60 tma (6.5) F is the two —dimensional Fourier transform of 0 over all 60 frequencies and a selected section of time. The variable u is the transform of frequency and the variable v is the transform of time. Since 0 is real, the spectrum of 0 is conjugate symmetric. Consequently, only half the complete two-dimensional transform need be computed in each direction to compute the power spectrum. The coherently averaged power spectrum, e(u, v), is the two-dimensional Fourier transform of the complete sub — set. In terms of F e(u, v) = IF(u, v, 0, 25611 2, u = 1,.., 30, v = 1f..., 128 (6.6) The incoherently averaged power spectrum, e' (u', v') is the average of the power spectra of the first, second,

119 third, and fourth time quarters of the 0 subset: E' (u, v ) 14 E F(u', v', 64k, 64 2 k=0 u' = 1,..., 30, v' = 1,..., 32 (6.7) The relation between the coherent (unprimed) and incoherent (primed) transform variables is u = u', u = 1,.., 30 (6.8) v = 4v', v' = 1,.., 32 (6.9) Thus the coherently and incoherently averaged spectra are comparable point for point in the u direction. In the v direction, the coherently averaged spectrum is computed on a net four times as fine as the incoherently averaged spectrum. The transform variables u and v are easily related to physical units. Consider first the coherently averaged spectrum. A peak at (u, v) implies there is significant amplitude present on a diagonal in the f-t plane which has slope 60 Af u- 1 At 256 v - 1

1'20 in f-t coordinates. Since 256 time points span 424.96 minutes and 60 frequency points span 50 Hz, the slope in physical units is v-1 50 v-l m 117.658 millihbertz/min. (6.10) c 424.96 u-l u-l Similarly, for the incoherently averaged spectrum, v'-l 50 v'-l m. = 1046 7 2470.632 ~ milliHertz/min. m 106.24 u'-l u'-l (6.11) Pictures of the coherently and incoherently averaged power spectra are shown in Fig. 6.17. Only the lower-frequency values along the v axis are shown, as the highfrequency values were essentially zero. The dc value, (u=l, v=l), in each analysis is suppressed to facilitate plotting. As expected, there are strong components on the v axis in both spectra. These arise from frequency-independent features in 0. There are also noticeable offaxis components, particularly in the area bounded by u' = 5,..., 7 and v' = 3,...,. A list of the locations of the largest of these off-axis peaks and their corresponding slopes is given in Table 6.1.

121 Coherently averaged V u V Incoherently averaged Fig..6.17. Two-dimensional power spectra of the phase rate 0 (f, t) from Fig. 6.15

122 Table 6.1. The location of off —axis peaks in the power spectrum of 0 and their corresponding slopes u' v' slope m. 5 5 0.470 Hz/min. 6 3 1.176 Hz/min. 3 3 1.647 Hz/min. These components correspond to the diagonal movements of 0 in the f-t plane mentioned in Section 6.4. 6.6 The Power Spectrum of Frequency-Independent 0 Features Frequency-independent features in e transform into spectral components lying on the v and v' axes in Fig. 6.17. Since these components are the strongest in the power spectra, they have been analyzed separately. Figure 6.18 is a graph of e(0, v), the coherently averaged spectrum along the v axis. As 256 points along the t axis were included in the transform, the period of any component is 256 -99.6 seconds v-1

123 212 mins 42.5 mins 24.28 mins 1000 J1 ~ dl17 mins, 15.73 mins 14.16 mins 100 7.86 mins 10 10 0 8 16 24 32 40 48 56 64 v-1 -- Fig. 6.18. Coherent power spectrum of the phase rate along the time axis

124 The periods of the more powerful components have been marked on Fig. 6.18. Of particular interest are the peaks at 17 and 15.73 minutes period, which are due to the 16-minute waves of Section 6.4. Also noticeable are peaks around 8 minutes period, the second harmonic of the 16-minute wave period. The presence of these second harmonics is strong evidence for the presence of the 16-minute waves in the original data. The graph of the incoherently averaged power spectrum, 8' (0, v'), in Fig. 6.19 shows many of the same features as Fig. 6.18. The graph is to the same scale as Fig. 6.18 along the v direction. Strong, harmonically related components are noticeable at v' = 7 and v' = 14, periods of 15.1 and 7.58 minutes respectively. Activity at lower frequencies is not nearly so prevalent in the incoherently averaged spectrum as in the coherently averaged spectrum. The shape of the graph in Fig. 6.19 suggests a model for the spectrum. The power as a function of frequency drops exponentially until a break frequency is reached around 10 minutes period. Above the break frequency, the power is approximately constant. This spectrum could well be due to the superposition of two physical processes: one with exponentially falling spectral density and one with constant, low-level spectral density. The coherently averaged power spectrum of Fig. 6.18

125 1000 100 1 0.i!,! J I U I I i, I 0 2 4 6 8 10 12 14 16 v'-1 ----- 7Fi?. 6.19. incoherent power spectrum of the phase rate along the time axis

126 was smoothed by grouping five adjacent frequency components to estimate the power at frequencies corresponding to the values of v' in the incoherently averaged power spectrum. Let va = 4-v'. Then the smoothed power estimate at v is e' a ea (0, va) - (0,va + k) (6.12) k=-2 where v runs from 1 to 128 in increments of 4. If the a argument va + k is undefined in Eq. 6.12, e is chosen to 1 be zero and the normalizing factor of - is altered accordingly. A plot of ea is shown in Fig. 6.2.0, again to the same frequency scale as Fig. 6.18. The first and second harmonic components due to the 16-minute waves are still present, although somewhat attenuated over Fig. 6.19. The general shape of the spectrum, however, is unaltered. 6.7 Conclusions The visual analyses of the phase of the wide bandpass transfer function phase have yielded the following insights into the propagation of underwater sound: The phase is not pathological. The simple process of assigning each phase measurement to the proper sheet of the Riemann Surface of the arctangent produces as highly structured phase surface.

127 1000 C 00 10 1.0 0.1 J ~ I | 0 8 16 24 32 40 43S 56 64 v -1 a Fic. 6.20. Five point movinc average of the coherent power spectrum of Fig. 6.18

128 T* he phase of the transfer function is coherent across frequency. The exceptions occur during propagation path structure changes and in regions of low signal-to-noise ratio. * Contours of constant phase move on diagonals in the f-t plane as predicted by Eq. 2.15. * Phase rate filtering is a useful method to counter-act the effects of low signal-to-noise ratio. * Noticeable diagonal movements of the phase rate are associated with low signal-to-noise ratio. * The presence in the phase rate of waves with 16minute periods is probably due to internal waves in the Straits of Florida. Further, these 16-minute waves are independent of frequency and signalto-noise ratio. Propagation path structure changes appear as large peaks in the phase rate, independent of signalto-noise ratio. The two-dimensional power spectrum of the phase rate confirms the movement of the phase rate along diagonals in the f-t plane. The power spectrum of the frequency-independent features in the phase rate falls exponentially to a period of 10 minutes. The power spectrum at frequencies from 10 minutes to 3 minutes period is approximately constant.

129 The following additional conclusions may be drawn regarding the anal-ysis techniques themselves: Any phase tracking algorithm which does not utilize signal-to-noise ratio information will malfunction at low signal-to-noise ratio. Amplitude information alone is not enough to rectify this problem. Noise levels must also be measured. Algorithmic techniques are preferable to heuristics to handle problems of features induced by noise. The combination of phase rate filtering with the ploy of forcing the phase rate to zero in low signal-to-noise ratio regions allows ready recognition of noise-induced features.

CHAPTER VII INTER-LINE CORRELATIONS The question is often asked, how far apart in frequency must two underwater sound receptions be so as to be essentially uncorrelated. Underlying this question is the intuition that the ocean is akin to a'first-order Markov process, where a single number, in this case the decorrelation bandwidth, can describe the correlation properties of receptions. It is shown below that the idea of decorrelation bandwidth cannot be applied to the wide bandpass receptions of this study, and further that the phase of receptions exercises almost peremptory control over their correlation. Given two complex-valued time series, x(t) and y(t), t = 1,... tmax define a time-series-valued function <<x, y>>(t) on x and y to compute the square of the normalized, sliding correlation coefficient: t+n | E x(i) y*(i) <<x, y>>(t) = i=t (7.1) t+n t+n x (i) x ~H I y(i) 2 i=t i=t t = 1,..., t + 1 - n 130

131 <<x, y>> is the square of the correlation coefficient between n terms of the two time series and ranges in value from 0 to 1. For this study, n was chosen to be 10. Thus 99.6-10 seconds = 16.6 minutes of data are correlated to produce one value in <<x, y>> Figure 7.1 shows p2 = <<S(f, t), S(31, t)>> the sliding correlation between all of S and the carrier line. If S(f, t) becomes decorrelated with the carrier S(31, t) with increasing distance from the carrier frequency, p2 should droop noticeably at the edges of the frequency band. However, there is no noticeable drooping of p2 at the band edges in Fig. 7.1. Indeed, no general dependency.of p2 on frequency manifests itself at all. Consequently, if a decorrelation bandwidth exists for this frequency band in the Straits of Florida, its value must be >>25 Hz. Further, Fig. 7.1 casts doubt on the existence of any decorrelation frequency within the bandwidth attainable by commonly available acoustic transducers. Most of the regions in p2 where the correlation is low can be recognized as regions of high path loss (cf., Fig. 5.3). Very apparent are the leftward drift in the front of the surface, the horseshoe feature, and the valley at the left rear. The dynamic range of the p2 surface is 35.5 dB. Since low signal-to-noise ratio appears to be strongly connected to low correlation values, the question may be

132 Fig. 7.1. The sliding correlation of S(31, t) with S(f, t)

133 asked how much of this dynamic range is due to raw amplitude changes and how much is due to phase uncertainties induced by low signal-to-noise ratio. To answer this question, sliding correlations were. done separately for the amplitude and for the phase of S(f, t). P2 = << I S(f, t) I, II S(31, t) I >>, the amplitude correlation, is shown in Fig. 7.2, and p~2 - e(frt) jO(31,t) p <<e(f et) the phase-only correlation, is shown in Fig. 7.3. The dynamic range of p2 is 7.06 dB a while that of p2 is 55.3 dB. Plainly, most of the dynap mic range and (from the pictures) a number of features of p2 are due to the phase measurements which were computed at low signal-to-noise ratio. In brief, noisy receptions have low correlation, while noise-free receptions are all highly correlated.

134 / i~~~~~/ Fig. 7.2. The sliding correlation of the amplitudes of S(31, t) and S(f, t)

135 Fig. 7.3. The sliding correlation of the phase portion of S(31, t) and S(f, t)

CHAPTER VIII CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK In the past, underwater sound propagation experiments utilizing acoustic transducer sound sources have been confined to measurements at one or two spot frequencies to reduce the data analysis problem to manageable proportions, despite the fact that multifrequency measurements held great promise for insight into propagation mechanisms. This study has produced data display tools which allow and indeed encourage simultaneous analysis of many frequencies. Using the techniques reported here, large amounts of multifrequencyqbandpass propagation data may be easily scanned by a human observer viewing computer-generated pictures. The artifice of defining transformations on the propagation data which turn regions of interest into peaks and edges in the pictures allows the observer to track these regions in both frequency and time, and thus to gain insights into propagation which are impossible to glean from single-frequency data. It is in the nature of studies such as this in which a powerful new tool is introduced to raise more questions than they answer. Nonetheless, the following areas seem to hold the most potential for future work: Generation of time-varying power spectrum, path loss, and phase rate pictures for Straits of Florida data taken over a much longer period than the 7.58 136

137 hours analyzed here. Generation of time-varying power spectrum, path loss, and phase rate pictures for data taken from other underwater sound test ranges, to look for similarities between their propagation and that of the Miami-Bimini range. Generation of time-varying power spectrum and path loss pictures for underwater ambient noise. Again, it cannot be iterated too strongly that in future work the noise level must be measured so that features caused by noise can be distinguished from features induced by the vagaries of propagation.

APPENDIX A THE SOFTWARE ENGINEE-RIJNG OF THE DISPLAY PROCEDURES Appendix A is concerned with the design of the display procedures from the view!point of modern software engineering techniques. The use of the display procedures themselves is presented in Appendices B and C. A.1. Implementation of the Display Programs The data surface display programs for parallel projection and contour mapping were each structured in a modular fashion [21],[22]. The kernel of each display program is a set of general purpose, device-independent display procedures. These procedures produce coordinate pairs rather than device commands for a digital plotter or an interactive display. This organization permits easy substitution of alternate plotting devices as they become available or especially suited to the investigation. Surrounding the display procedures are routines which ~ communicate with the user concerning picture oarameters, ~ read the data surface coordinates into memory, ~ invoke the display procedures, e derive plotting device commands from picture vectors, and produce labels on the pictures. 138

139 A.2. Data Representation The S(f, t) data points were taken on a regular grid: equal steps in time and equal steps in frequency. The internal computer representation of S could consist simply of the array of z values and two spatial increment values, one for the x axis and the other for the y axis. However, recognizing that the x- and y-axis data locations need not be uniformly spaced for all potential applications of the display procedures leads to a more general, parametric data representation. In this representation, actual x and y coordinate values are stored in individual vectors such that if Z = S(f0, t0), then inside the display procedures, Z(I, J) = S(X(I), Y(J)) (A.1) X(I) f (A.2) Y(J) = to (A.3) This combination of one-dimensional X and Y arrays and a two-dimensional Z array is then the machine representation of the data surface to be displayed. A.3 Development Goals The goals of the display procedure development process were to produce procedures which were Verifiably correct.

140 Fast in execution, Flexible and easy to use, and * Capable of sequential realization for on —line use. These are discussed in more detail below. In addition to the above goals, the author desired to build on previous work in the field, and thus the new procedures were to be refinements insofar as possible of extant procedures. The parallel projection procedures were built upon a set of subroutines available in the library of the University of British Columbia Computing Centre [16]. The contour mapping procedures were built upon an algorithm of Hartwig [17] and information in a paper by Walters [14]. A.3.1. Verification of Correctness. The wisdom of formally arguing for the correctness of a program has recently been recognized ([18] [19] [20]). Even if a formal argument is never made, thinking about such arguments focuses the implementor's mind on the error-prone input/ output relationships of basic blocks in his program and thus offers a path to the production of high-quality programs. Since correctness arguments are facilitated by the use of high-level programming languages, it was decided that only such languages would be used to write the display procedures. A.3.2. Speed of Execution. During the initial development of a procedure, speed of execution should be

141 secondary to considerations of correctness, readability, and maintainability. Too many implementors interchange these priorities, not realizing that efficiency can only be a property of a CORRECT program. Experience emphatically shows that during the development of a program the most substantial gains in execution speed are attained through changes of algorithm rather than through exploitation of specific properties of a particular machine. Such changes are also facilitated by the use of a high-level language. Given these considerations, the development plan when writing the display procedures was as follows: First, write correct procedures in a high-level language and argue their correctness. Next, using measurements of the central processor utilization of each program statement, replace those algorithms which prove computationally inefficient by faster algorithms. This replacement process iterates until speed goals are met. Provided the input/output properties of the faster algorithms are the same as those of the slower algorithms they replace, the global correctness proof of the procedure remains unchanged. Only the new algorithm need be proved. A.3.3. Flexibility. By flexibility is meant yielding to the user complete control over the behavior of the display procedures. There should be no significant aspect of the output picture or display procedure behavior which cannot be controlled by a parameter in the calling sequence.

142 One way to achieve this goal, recently presented in the literature, is to make the set of display procedures be a "heap" (Spier, [23]). A heap is a memoryless collection of subsystems in which all actions taken by the collection are controlled solely by the variables input to the collection. To make a set of procedures be a heap, the parameters which determine the course of the computation (i.e., the observable plant-input state variables) must be placed in the calling sequence of the procedures. A.3.4. Sequential Realization. Given their success in off-line analyses, the display procedures will certainly be useful in the on-site, real time data acquisition environment. In this environment, the data surfaces are not available en mass, but rather are gathered piecemeal by the real time data acquisition system. Display generation algorithms would be best matched to such an environment if they were capable of sequential realization, i.e., if they also could accept their input data piecemeal and generate pictures incrementally. Feedback to the experimenter from displays could thus be provided in real time. A.4 Achievement of Design Goals in the Parallel Projection Procedures The correctness performance of the parallel projection procedures has been encouraging, although not perfect. Correctness was argued extensively while the procedures were being written, and contributed significantly to the

143 overall understanding of the algorithm and of the problem it was written to solve. Nonetheless, the procedures do have a minor malfunction which is due to an error in the specification (i.e., the procedures implement their specification perfectly insofar as is known). Fortuitously, the malfunction manifests itself in an obvious way to the user so that he knows that the output is incorrect. Nonetheless, such a mistake compromises the trustworthiness of the procedures to the end user and cannot be lightly disregarded. The effort to prove the correctness of the parallel projection procedures was worthwhile and produced a number of insights: The process of applying top-down design (Dijkstra, [20 ) with argument for correctness at each level should extend to the problem specification as well as to the code for the algorithm. Confirming Dijkstra's contention, borrowing another programmer's coded algorithm can be detrimental to correctness, especially if correctness is not argued for the routine being borrowed. Although the changes to be made may seem minor, and although the algorithm may seem to be quite good, many pitfalls will be encountered in the process of adapting the other's algorithm to one's own situation. One must actively husband correctness when building on previously written code.

144 The process of deriving the correctness assertions for an algorithm improves the programmer's comprehension of the algorithm, making a correct, computationally efficient implementation more probable. * Inclusion of the correctness assertions as comments in the body of the code enhances readability and maintainability. * Some method of separating the correctness assertions from more mundane comments in the code would improve program comprehensibility. * The choice of programming language impacts substantially the effort needed to prove correctness. For example, the author feels that the use of PASCAL [24] [25] would have cut the effort involved in correctness arguments by a substantial margin. The goal of sufficient speed in execution has been met. For example, the algorithm will produce coordinates for a view of a 16,000-point data surface (S(f, t)) in 24 cpu seconds on an IBM/360 —67 computer running under the MTS operating system. The hidden line elimination algorithm consumes 70 per cent of this time. Sixteen additional seconds are required to read the data, produce digital plotter commands from the coordinate specifications, and do other miscellaneous operations. The hidden line elimination algorithm of Williamson [13] would undoubtedly cut the total execution time appreciably, but it was unknown to the author at the time the procedures were written.

145 Nonetheless, the parallel projection procedures are sufficiently inexpensive to use that the author has not cut an investigation short on economic grounds. In addition, other users who have switched to the author's procedures from a procedure which ran eight times slower have remarked that they feel a freedom to investigate their data which was not present under the older procedure. The parallel projection procedures were written before the author became familiar with the heap concept. Thus all the state variables which should be present in the calling sequence are not there. In particular, assumptions are made about error handling over which the user should have control. An additional problem arose, involved with making the total set of input state variables available to the user through the calling sequence. One consideration is that the user should not have to specify (and the machine should not have to handle) those input variables which will change little or not at all between calls. Conversely, Spier [23] can be interpreted to ask for passing the total input state vector to the procedure in one call. These two points conflict on a number of grounds. This author has opted to give the parallel projection procedures a tree — like hierarchy of entry points; the topmost entry points set those state variables which change slowly while the bottommost set those state variables which change fastest. Checking is incorporated to ensure that no procedure lower

146 in the hierarchy is allowed to run unless its predecessor higher in the hierarchy has already run successfully. The various entry points can then be viewed as a set of interlocked processes which could be run in parallel. It is significant for the design of future computer programming languages that an operating system issue appears in this, a strictly applications-oriented procedure. Although the current implementation of the parallel projection procedures requires that all the points in the data surface be in memory simultaneously, a sequential implementation is possible. (The hidden line eliminator already operates sequentially.) Once the user has specified rotation and elevation angles and scaling parameters, processing may proceed sequentially, handling one vector of data samples along either the x or y axis at a time. A.5 Achievement of Design Goals in the Contour Mapping Procedures As with the parallel projection procedures, the correctness performance of the contour mapping procedures has been quite encouraging albeit not perfect. Hartwig's procedure contains a "don't care" condition because of the interface he used to the digital plotter. In attempting to generalize his procedure and make it device independent, the author glided over the "don't care" and thereby introduced an error. Had Hartwig's original code been in a language which had a CASE statement, the error would likely

147 not have occurred —it derived from not handling all members of a dense subset of the positive integers. The programming language PASCAL would have particularly been at home handling this situation: the preferred style in the language would never have even let the "don't care" arise. Hartwig's algorithm is at once both more easily understood and yet more computationally complex than the parallel projection algorithm. Consequently, the author did not include correctness assertions as comments to the code, relying instead on other explanations. of the algorithm for documentation of the intent of the code. The contour mapping procedures meet the goal of fast execution. The algorithm will produce coordinates for the contour map of a 16,000-point data surface, contoured at eleven levels, in 12 cpu seconds on an IBM/360-67 computer running under the MlTS operating system. A number of heuristics are applied to reach this speed; based upon measurement, Hartwig's procedure takes one-third more time to execute than the procedures developed here. The total time, including overhead to produce a contour map, is just about that of producing a parallel projection —40 cpu seconds. Again, the author felt no economic pressure to terminate his contouring investigations. The contouring procedures are memoryless, and therefore must be a heap in the sense of Spier. Thus the user has access to all the input state-controlling variables

148 in the procedures, and the goal of flexibility has been achieved. The current implementation of the contour mapping procedures is essentially sequential. The implication in Appendix C that much of the data to be contoured should be in memory at one time has been inserted to discourage repeated, high-overhead calls to the procedures. In fact, the contribution of Hartwig's algorithm (about which, unfortunately, he says nothing) is that a sequential realization of a contouring algorithm seems invariably to possess lower computational complexity than algorithms which follow contours around on the data surface. A.6. Conclusions Correctness, flexibility, modularity, and speed were the main implementation goals for the display procedures. Correctness was ensured by program verification in the sense of Hoare [18]. Flexibility and modularity were achieved through application of the software engineering principles enunciated by Dijkstra [20], Parnas [21], and Spier [23]. Speed was attained through the selective replacement of slow algorithms by faster algorithms. The final products are demonstrably correct simple to use, and acceptably fast.

APPENDIX B DISPLAY PROCEDURES FOR PARALLEL PROJECTION VIEWS OF DATA SURFACES Type of Routine. A set of Fortran IV SUBROUTINE subprograms, collectively called PERSYS. The public entry points are PERS, PERSP, and PPSET. Purpose: The PERSYS subroutines will produce coordinates for a parallel projection view of a surface which is defined on a regular grid. Portions of the data surface which are hidden from the viewer are not drawn. Surfaces which are represented hby scattered data points or non-rectangular arids may not be drawn directly by these subroutines. How to Use: PPSET is called once before any calls to either PERS or PERSP to initialize the routines; once PPSET has been called, it need not be called again unless the user wants to change the value of one (or more) of the PPSET param — eters. PERS is called to produce picture coordinates for data surfaces defined on regular grids. Since PERS will be the most commonly used entry point, its use will be covered here. The use of PERSP is covered at the end of this Appendix. (PERS is in actuality just a preprocessor for PERSP.) 149

150 No matter whether PERS or PERSP is used to produce coordinates, PPSET must be called first to initialize the system. The program then reads like the following. REAL Z(IDX,?)X(IX), Y(IY) EXTERNAL PENUPS, PEIUDNS LOGICAL SWCHES(5) CALL PPSET(PENUPS, PENDNS, IRESOL, HSIZE, VSIZE) CALL PERS(XY, IX, Y, IY, Z, IDX, ABOUT, ABOVE, RY2X, RZ2X, SWCHES) The parameters to the PPSET call are: PENUPS —The name of a subroutine which PERSYS can call to dispose of an (x,y) picture coordinate pair. The subroutine should cause the beam or pen to be positioned invisibly to the point (x,y). x and y:are of type REAL. PENDNS —The name of the subroutine which PERSYS can call to dispose of an (x,y) picture coordinate pair. The subroutine should cause a visible line to be drawn from the beam or pen position prior to the call to the point (x,y). x and y are of type REAL. IRESOL —An integer constant or an INTEGER*4 variable which specifies the horizontal resolution of the hidden line eliminator in PERSYS. The value of IRESOL is tied closely to the device upon which the picture will be displayed. For example, on a display terminal with 1024 displayable horizontal points, an

151 IRESOL value of 1024 would cause at most a onepoint error in the horizontal position of a line; specifying IRESOL to be 2048 would remove even this one-point error, at the expense of doubled processing time. If the output device is a plotter, IRESOL should be specified such that the resolution of the hidden line eliminator is on the order of the horizontal resolution of the plotter. If the picture to be drawn is spatially coarse (i.e., has no lines of very short length), the user may specify lower resolution values in order to speed hidden line processing (the cost of which is directly proportional to IRESOL); this may, however, introduce an unacceptable level of distortion into the picture. (A vertical resolution of 32768 is provided automatically and at no added cost by the hidden line eliminator.) HSIZE, VSIZE —Two REAL numbers which fix the scaling of the output coordinates. The picture coordinates will be scaled such that no output point lies outside the rectangle bounded by the points (0., 0.) and (HSIZE, VSIZE) The parameters to the PERS call are: X A work vector of type REAL dimensioned at least IX entries long in the calling program. IX An integer constant or an INTEGER*4 variable giving

152 the number of points along the first dimension of the Z array which are to be included in the picture. IX must be > 2. Y A work vector of type REAL dimensioned at least IY locations long in the calling program. IY An integer constant or an INTEGER*4 variable giving the number of points along the second dimension of the Z array to be included in the picture. IY must be > 2. 7Z A doubly-dimensioned array of real numbers which contains the z coordinate values representing the surface. Each number in array Z gives a height associated with one position in the x-y plane, as discussed in Appendix A. Z(I,J) gives a value for the surface at the Ith grid line in the xdirection and the Jth grid line in the y-direction. The actual shape of the grid is given by argument RY2X (described below). NOTE: The array Z is destroyed by the subroutines. IDX An integer constant or an INTEGER*4 variable which must be equal to the value of the first subscript of the array Z as set in the DIMENSION statement. If Z is declared: DIMENSION Z(52,30) then: IDX = 52 The grid itself need not be the same size as array

153 Z although it may not be larger (i.e., IDX must be > IX). The actual size of the grid is given by the two arguments IX and IY. The restriction IDX > 2 is imposed by PERSYS for error-checking purposes. NOTE. An incorrect IDX is a major source of error. IDX need not, in general, be equal to the value of IX. If the plot contains unusual cliffs which don't belong, the cause is usually an incorrect IDX. (The viewing angle or position of the observer is determined by the next two arguments.) ABOUT A real number or a REAL*4 variable which gives the viewing angle in degrees measured counter-clockwise (to the right) from the negative y-z plane. ABOUT may be any angle between 0.0 and 89.0 degrees. If it is desired to view the surface from some other angle, it is necessary to re-enter the data into array Z in a different orientation. ABOVE A real constant or a REAL*4 variable which gives the viewing angle in degrees above the x-y plane. ABOVE may be any angle between 0.0 and 89.0 degrees. (The next two arguments describe the actual physical shape and size of the grid superimposed on the surface.) RY2X A real constant or a REAL*4 variable which gives the ratio of the length of the y-side of a grid

154 rectangle to the length of the x-side. Consequently, if the grid is square, RY2Z = 1.0. If the grid is rectangular with the y-side of the rectangle only one half of the length of the x - side, then RY2Z = 0.5/1.0 = 0.5. If the grid is not uniform for all data points, consult the information on entry point PERSP. To keep the user from doing something grossly in error, an error check' to ensure that 0.02 < RY2X < 50 is made inside PERSYS. RZ2X A real constant or a REAL*4 variable which gives the ratio of the maximum value of the surface minus the minimum value of the surface to the length of the surface in the x-direction. This variable determines whether the resulting plot will look like a low rolling plane (RZ2X small) or a range of mountains (RZ2X large). Since the purpose of this program is primarily to show the shape of the surface, and since height measurements cannot be taken from the plot, this variable may often be approximate or be selected to produce pleasing results. To keep the user from doing something grossly in error, an error check to ensure that 0.02 < RZ2X < 50 is made inside PERSYS. SWCHES A five-element LOGICAL vector used to control various picture options. The assignment of values is

5bb5 as follows:' SWCHES(1) T: Baselines are to be drawn. F: No baselines are to be drawn. SWCHES(2) (Valid only if SWCHES(1) is true). T: Baseline height is at Z = 0 F: Baseline height is at the minimum Z in the data surface. SWCHES(3) T: Write an informational message about various PERSYS parameters (IX,IY,IDX,ABOUT,ABOVE, HSIZE,VSIZE) on logical I/O unit 10. F: Write nothing on logical I/O unit 10. SWC}T{:E (4) Tz Draw a line in the picture between Z(IJ) and Z(I,J+l) for J = 1...,IY=l,... IX. F. Do not draiw the immediately above —mentioned lines. SWCHES(5) T: Draw a line in the picture between Z(IJ) and Z(I,J+l) for J = 1a...,IY=1 and I=1 and I=IX only, F: Do not draw the immediately above-mentioned lines. Method: PERSYS consists of three subroutines. PERS calls subroutine PERSP (described further on) which in turn calls HIDLIN which does hidden line elimination and prod — uces coordinate pairs. Subroutine PERS sets up arrays X and Y required by PERSP and scales the data according to the argumentRZ2X. Subroutine PERSP rotates the data surface about the z and x axes and forms the x —z parallel projection of the rotated surface; array Z is used to store intermediate data to reduce the number of computations re — quired. The projection is then modified by assigning each point in the projection to the nearest intersection of a

156 discrete grid IRESOL points wide by 327,68/;points high. If RIRESOL "s sufficiently large, this modification will dis — tort the picture insignificantly; hidden, line elimination, however, is greatly facilitated by this modification. PERSP plots'baselines if requested, and, then plots the modified projecton through subroutine HIDLIN which contains the coding to suppress the plotting of hidden lines. HIDLIN dynamically allocates an array'according t6''the ]valuelDof IRESOL''s part of the hidden line processing. Logical TO/Units,Referenced: 3 Error comments from PERSYS if parameters are invalid. 10 Informational output from PERSYS if SWCHES(3) is set.to true. Memory Requirements: PERS 1664 by..tes PERSP~ -'4498 bytes HIDLIN 2284 byte-s~ 2*IRESOL bytes for hidden line elimination Error Handling: I f any of- the parameters IX, IY, IDX, RY2X, RZ2X, ABOUT,-'iand'ABOVE is outside its valid range, PERSYS will write an'error comment to logical I/O unit' 8, the picture wil'l not be produced, and'an immediate return to the user wille be madea. If the user entered via. PERS, a simple RETURN is made when an error is encountered, while if the

157 user entered via PERSP, a RETURN 1 is executed when an error is encountered. Subroutine PERSP Purpose: Subroutine PERSP may be used to plot a perspective view of a surface when the data grid is not necessarily composed of rectangles all of the same size but where the grid is still laid out in a rectangular manner. T.o X7. to Uses This subroutine is called as- a normal FORTRAN subroutine subprogram. CALL PERSP (XP,IX,YP,IY, IDX,ABOUT,ABOVE,SWICHES) IX,IY,Z,IDX,ABOU2,ABOVE, and SWCHES are the same variables as described in the writeup for subroutine PERS. Since the grid in this case is not necessarily regular, the arguments RY2X and RZ2X have been deleted and the shape and size of the grid is described by the arrays XP and YP. XP is a REAL vector, dimensioned at least IX in length, which contains the x coordinates of the grid in the same scale as the data in the Z array. YP is a REAL vector, dimensioned at least IY in length, which contains the same information for the y coordinates as XP specifies for the x coordinates.

158 The grid to be plotted need not be regular. The grid corners are located by the values of the XP and YP arrays as discussed in Appendix A.

APPEND IX C CONTOUR MAPPING PROCEDURES FOR DATA SURFACES Type of Routine: CONSYS is thle collective name for the contour mapping procedures, which are a collection of IBM Fortran IV SUBROUTINE and FUNCTION subprograms. The public entry points are CONTUR, CONSET., CONLBL, ATND CONEND. The routines are not in ANSI standard Fortran IV, and thus would need some work before they could be transported to a non-IBM installation. Moreover, the routines do dynamic storage allocation of some work spaces and thus are somewhat MTS-dependent; this latter objection to transportability could be removed by'fixing the dimensions of the work arrays and recompiling the source modules. Facilities Available: The CONSYS subroutines will produce coordinates for a contour map of a data surface which is defined on a rec — tangular grid. Surfaces which are represented by scattered data points or nonrectangular grids may not be directly contoured by these routines. CONSYS contains a primitive contour labeling facility which the user may invoke by setting a switch in the subrouting calling parameters. For the applications envisaged for CONSYS (large data surfaces), the cpu time necessary to produce "nice' labels was considered an extravagance. 159

160 (The CPU time would increase by a factor of at least 1.5.) CONSYS does not require that a user have the data to be contoured all present in memory at one time. Thus the user may produce very large maps by calling CONTUR repeatedly, once for each different region of the data surface. The user has control over the drawing of border lines for the map and can thus omit the border line when another map section is to be adjacent to the one currently being produced. Since CONSYS dynamically allocates work areas, there is effectively no limit to the number of data points in the data surface to be contoured. To keep CONSYS device independent, the user must write a labeling routine if he/she desires to invoke the labeling option of CONSYS. This routine is often only 10-20 statements long, but is usually dependent upon the graphical output device to which the user sends the map coordinates. How to Use: The user must call subroutine CONSET once to initialize CONSYS before calling upon CONTUR to produce contour line coordinates. Once CONSET has been called, it need not be called again unless the user wants to change the value of one (or more) of the CORSET parameters. If the user desires to use the labeling facility in CONSYS, he/she must call subroutine CONLBL to pass in various labeling parameters. CONLBL must be called after

161 CONSET has been called, but before the first call to CONTUR which specifies that labeling is to be done. The user may then call on subroutine CONTUR to produce coordinates for a contour map. The coordinates produced are in the- coordinate system and value range which the user specifies in the parameters to his/her call to CONTUR; thus the user may make repeated calls to CONTUR either to produce new contours on an existing map or to map new data regions adjacent to those already mapped. When the user is finished producing contours, he/she must call the routine CONEND to release dynamically acquired storage. Calling CONaET: None of the calling parameters to CON4SET are changed by CONSET in any fashion. The calling sequence for CONSET is CALL CONSET(PNUP, PNDN, IOUNIT) where: PNUP is the name of a subroutine which CONSYS can call to dispose of an (x,y) contour line coordinate pair. The subroutine should cause the beam or pen to be positioned invisibly to the point (x,y) in the user's coordinate space. x and y are of type REAL. PNUP must be declared EXTERNAL in the calling program.

162 PNDN is the name of a subroutine which CONSYS can call to dispose of an (x,y) contour line coordinate pair. The subroutine should cause a visible line to be drawn from the beam or pen position prior to the call to the point (x,y) in the user's coordinate spce. x and y are of type REAL. PNDN must be declared EXTERNAL in the calling program. IOUNIT is an integer constant or an INTEGER*4 expresion which specifies the Fortran data set reference number (DSRLC.??) on which error comments from CONSYS are to be written. If the value of IOUNIT is less than zero or greater than nineteen (currently the largest legal DSRUN in MTS Fortran), then error comments will not be written. Conversely, if IOUNTI' is a legal MTS Fortran DSRN, then in the event of errors, error comments will be produced via FortranWRITE statements which look like WRITE(IOUNIT, format).. Contour Labeling: CONSYS has a primitive contour labeling facility which may prove adequate for a large number of applications. Wdith this facility, the user specifies in his coordinate system a number of lines of constant x and a number of lines of

163 constant y. Whenever a contour enters a grid cell which contains one of these constant-coordinate lines, the contour becomes a candidate for labeling in that cell. If the center of the grid cell is more than a user-specified distance away from the the centers of all other grid cells which have already been labeled, the contour is declared a successful candidate for labeling. In this event, CONSYS then calls a user-specified subroutine with three parameters: the center coordinates of the grid cell (x and y), and the z-value of the contour which entered the grid' cell; the user-specified routine is then responsible for producing a label for the contour and sending the label to the graphical output device. The values of the constant coordinate lines, the name of the labeling routine, and the distance parameter for declaring a labeling candidate to be successful are passed to CONSYS by calling CONLBL. Calling CONLBL: If the user does not wish to use the labeling facility of CONSYS, then he/she does not need to call CONLBL or to read this section. None of the calling parameters to CONLBL is altered in any fashion by CONLBL. The calling sequence for CONLBL is CALL CONLBL(LABELR, DIST, NXL, XLOC, NYL, YLOC) where LABELR is the name of a subroutine which CONSYS will

164 call to produce a label for a contour line. LABELR will be called via CALL LABELR(XC, YC, Z) where XC and YC are the type REAL x- and y-coordinates of the center of the grid cell containing the contour to be labeled. Z is the type REAL value of the z-coordinate of the contour line which is to be labeled. LABELR must be declared EXTERNAL in the calling program. DIST is a REAL constant or a REAL expression giving'the minimum distance in the user' s coordinate space upon which declaring a labeling candidate successful will be conditioned. If the distance from the center of a cell which is a candidate for labeling to any other cell which has already been labeled is greater than DIST units, then the candidate cell will be labeled by calling LABELR. If each label takes at most N characters with width S coordinate units per character, then a good value for DIST is S (N+1). DIST must be > 0. NXL is an integer constant or an INTEGER*4 expression giving the number of constant-x lines to be used in the labeling process. NXL must be > 0.

165 XLOC is a REAL vector at least NXL entries long. (An argument must be supplied even if NtXL is zero.) Each entry in XLOC is the x coordinate in the user's coordinate space of a line of constant x value (i.e., a vertical line) to be used during the labeling process to determine candidate cells for labeling. (See "Contour Labeling:; above.) NYL is an integer constant or an INTEGER*4 expression giving the number of constant-y lines to be used in the labeling process. NYL must be > 0. YLOC is a REAL vector at least NYL entries long. (An argument must be supplied even if NYL is zero.) Each entry in YLOC is the y coordinate in the user's coordinate space of a line of constant y value (i.e., a horizontal line) to be used during the labeling process to determine candidate cells for labeling. (See "Contour Labeling" above.) Note: NXL and NYL may not both be zero at the same time. Calling CONTUR: Once COUNSYS has been initialized by a call to CORJSET (and, optionally, a call to CONLBL), the user may call CONTUR to produce coordinates for contour lines. Jone of the calling parameters to CONTUR are changed in any way by

166 CONTUR. The calling sequence for COIJTUR is CALL CONTUR(X, IX, Y, IY, Z, IDX, CZ, NC, SWCHES) where X is a vector of x coordinates at least IX entries long. Each coordinate marks a vertical line which is the vertical boundary of a grid cell. The coordinates in X are in the user's coordinate space. IX is an integer constant or an INTEGER*4 expression giving the number of x coordinate values in the X vector to use, i.e., the number of points along the first dimension of the Z array to be considered in the contouring process. IX must be in the range 2 < IX < IDX. Y is a vector of y coordinates at least IY entries long. Each coordinate marks a horizontal line which forms the horizontal boundary of a grid cell. The coordinates in Y are in the user's coordinate space. IY is an integer constant or an INTEGER*4 expression giving the number of y coordinate values in the Y vector to use, i.e., the number of points along the second dimension of the Z array to be considered in the contouring process. IY must be > 2.

167 Note: The data grid upon which contours are to be produced need not be regular. The corners of the grid are located by the values in the X and Y arrays, as indicated in Appendix A. Z is a two-dimensional array of REAL values which are the z coordinates representing the data surface. The first dimension of Z is IDX while the second dimension of Z must be > IY. Each number in array Z gives a height associated with one position in the x-y plane. Z(I,J) gives a value for the surface at the Ith grid in the x — direction and the Jth grid line in the y-direction; i.e., the correspondence between Z and the X and Y vectors is that Z(I,J) = f(X(I), Y(J)) IDX an integer constant or an INTEGER*4 expression which must be equal to the value of the first subscript of the array Z as set in the DIIMENSION statement. If Z is declared: DIMENSION Z(52, 30) then: IDX = 52 The data grid itself need not be the same size as array Z although it may not be larger (i.e., IDvX must be > IX). Thle actual size of the grid is given by the trwn arguments IX and IY. The restriction IDX > 2 is imposed by CONSYS for

168 error-checking purposes. Note: An incorrect IDX is a major source of error. IDX need not, in general, be equal to the value of IX. If the map contains unusual jagged contours which don't belong, the cause is usually an incorrect IDX. CZ is a vector of REAL values, at least NC entries long, which contains the z-coordinate values at which contours are to be produced. The coordinates in CZ are in the user's coordinate space. NC is an integer constant or an INTEGER*4 expression whose value is the number of values in the CZ vector for which contours are to be produced. NC must be > 1; although CONSYS will work for NC = 1, it will function more efficiently on a per-contour basis when NC is greater than one. SWCHES is a type LOGICAL vector having five elements used to control various contouring options. The assignment of values is as follows: SWCHES(1) T: Draw a border line between (X(1), Y(1), and (X(IX), Y(1)). F: Do not draw the above border line. SWCHES(2) T: Draw a border line between (X(1), Y(1), and (X(1), Y(IY)). F: Do not draw the above border line. SWCHES(3) T: Draw a border line between (X(1), Y(IY)) and (X(IX), Y(IY))

169 F: Do not draw the above border line. SWCHES(4) T: Draw a border line between (X(IX), Y(IY)) and (X(IX), Y(1) ). F: Do not draw the above border line. SWCHES(5) T: Label contours using the labeling scheme described above. F: Do not label contours. Calling CONEND: CONEND is called when all contouring of a specific surface has been completed. It deallocates the dynamically acquired arrays which were used in the labeling process during calls to CONTUR to remember the center coordinates of cells which were labeled. (These arrays are not automatically deallocated when CONTUR returns, thus insuring that successive calls to CONlTUR will not place labels too closely together.) If CONSYS was used with labeling turned off (i.e. SWCHES(5) is.FALSE.), then CONEND need not be called; if it is called, the call will be ignored and no diagnostic will be generated. The parameters to CONEND are used to return information to the caller in three variables regarding the perforance of the dynamic storage routines for labeling. The calling sequence is CALL CONEND (NREGEN, H-IAXSIZ, NUSED) where

170 NREGEN is the INTEGER number of regenerations of dynamic storage which took place. MAXSIZ is the INTEGER number of bytes of storage allocated by the labeling routines NUSED is the INTEGER number of bytes of storage actually used by the labeling routines. All three CONEND parameters must be integer variable names; they may not be integer constants or integer expressions. The user may choose to do nothing with the values returned in these variables, but if he calls CONEND, they must be supplied. Logical I/O Units Referenced: CONSYS references the I/O unit (Data Set Reference Number) IOUNIT passed into CONSET. If CONSYS must produce an error comment and CONSET has not been called (this is caused by calling the routines out of order), the unit 19 is used. Memory Requirements: CONTUR 4644 bytes DISTOK 1494 bytes CTQQ 2368 bytes CONCOM 56 bytes 2*NC bytes for sorting contour values 8 bytes per label Error Handling: Extensive checks are made to see if parameters supplied by the user are within prescribed ranges. If an

171 error is detected, then if IOUNIT is within the legal MTS range for logical I/O units, an error comment is written on DSRN "IOUNIT". A nonzero RETURN statement is then exe — cuted and the call is ignored. A summary of the error handling is given below. Module Type of Return Cause CONLBL RETURN 1 CONSET not yet called. RETURN 2 DIST, NXL, or NYL out of range. CONTUR RETURN 1 CONSET not yet called. RETURN 2 SWCHES (5) set, but CONLBL not called. RETURN 3 IX, IY, IDX, or 1NC out of range. Method: CONSYS uses an algorithm which is a modification of the algorithm of Hartwig [17]. The modifications appear to have cut the time required to contour large surfaces by a factor of two. The modified algorithm reads approximately as follows:

172 Sort the values in CZ into ascending order; for each cell in the grid do; for each contour value in CZ do; if the cell intersects the contour value at all then do; using the fact that CZ is sorted, find the smallest and largest values in CZ which do intersect the cell; call CTQQ (which utilizes Hartwig's algorithm) to produce the contour coordinates for the cell. end: end; end; The speedup over Hartwig's original algorithm is attained in the step "if the cell intersects the contour value at all," wherein the cell is very quickly rejected if no contours pass through it; this is the case most of the time. In the subroutine CTQQ, contour coordinates are produced for a rectangular cell according to the Hartwig method. Briefly, the center coordinates of the cell are computed as averages of the corner coordinates. The corners and the center are then connected with line segments which are in turn numbered 1 through 8 in a clock-wise direction:

173 3 (x,y ) 7 2x'yz, i1,j, l Each segment is then tested in turn to see if the required contour intersects with it. If there is an intersection,then linear interpolation along the line segment is used to fined the x and y coordinates of the intersection, and the coordinates are saved temporarily in an array. After all eight line segments have been processed in this way, the coordinate values which were saved are sorted into proper orcder to form contour lines and are output ]ny calling the user specified output routines. Below is an

174 an example of a contour line crossing line segments 1, 2, 4 and 5.

APPENDIX D BIRDSALL'S PHASE RECONSTRUCTION ALGORITHM The algorithm used to reconstruct the phase of S(f, t) was developed by T. G. Birdsall. It reconstructs the phase at frequency f independently of the behavior of S(f, t) at any other frequency. O(f, t) from Eq. 6.1 was chosen to lie in the range from 0 to 1 rev. Let 0(f, t) be the reconstructed phase angle whose value ranges over the real line. Since c(f, t) is always a fraction, reconstruction simply involves choosing an appropriate integer part I(f, t) for e(f, t) = I(f, t) + ~(f, t) (D.1) Let E 1 and 2 be defined by 1l = le(f,t-1) - e(f, t)| (D.2) E2 = le(f,t-2) - 6(f, t)| The reconstruction algorithm chooses I(f, t) to minimize the larger of C1 and ~2 2 The algorithm uses no information from adjacent frequencies when constructing the phase for frequency f, in spite of the relation between 175

176 frequencies expressed by Eq. 2.15. The computational form of the reconstruction algorithm is shown in Figure D.1. For initial conditions, G(f, -1) = e(f, 0) = 0 are used.

177 {This algorithm is presented in the } {programming language PASCAL [24] [25]} const fmax = 61; tmax = 274; type frequencyindex = 1..fmax; timeindex = l..tmax; var a, c, m: real; f: frequencyindex; t: timeindex; 6,4: array [frequencyindex, timeindex] of real; begin {trunc(x) returns the integer part of x} a:= f, t] + trunc(9[f, t-l]); if e[f, t-l] < 0. then a:= a - 1.; m:= e[f, t-2] - a + e[f, t-l] - a; if abs(m) < 1. then c:= 0. else if m < -1. then c:= 1. else c:=-1.; 6[f, t]:= a + c; end Fig. D.1. Birdsall's phase angle reconstruction algorithm

APPENDIX E THE PROBABILITY DENSITY ON PHASE AS A FUNCTION 01i' SIGNAL TO NOISE RATIO In signal space, assume a reception consists of a signal Sejf in two dimensional additive Gaussian noise having variance 2 in both the x and y directions. The probability density function of the reception in Cartesian coordinates is -[(x-S cos p)2]+(y-S sin p)2]! 202 (E.1) f(x,y) = e 2Tac2 Without loss of generality, we may let q = 0. Then -[ (x-S) 2+y2 ] 1 2a2 f(x,y) = e [x2+y 22xS+S2] (E2) 2782 2Tr 2 178

179 Let v2 be the signal to noise ratio in the plane: 2 S2 S 2 = S v = (E.3) 2 2 Then _ [x2+y 2_ /Zxv + v2] f(x,y) e(E.4) 2jr2 Changing to polar coordinates (r,e) with Jacobian r yields r2 /2vr cos 2] f (r1) - r e L22 G (E.5) f(r~2) e Completing the square in the exponent gives f reV2 [22 - cos + v2cos2]+v2cos2 -V 202 2r r e~2+v2cos26 re v cos] (E.6) 2C 2 27T

180 This density may be integrated over r from zero to infinity to yield f(e): -v - 1 +R-v/cosel f(s) e re dr (E.7) 0 Changing variables in the integral, let u = r - /2v cose (E.8) Then 00 2 -v +v2cos20 - f(e) = (cu+/2vcos6) e du (E.9) -/2vcose Breaking the integral into the sum of two integrals and performing the integration yields v 2+V2COS20 1 -v2cos 2 / 2 vcos6 f(0) = e v c- e + 2 D (Y/-2VCOS 8) (E.10)

181 where x U2 1 2 D (x) e du (E. 11) - 00 is the probability distribution function of the Gaussian density with zero mean and unit variance. Factoring and combining terms yields f(s) = e [1 + /D2v cos (2v ose (E.12) 2Tr -v2 cos2e Recognizing that -v2Cos2 is the reciprocal of the Gaussian density function N(x) with zero mean and unit variance with x = vv cose allows a further simplification: f(s) = 2 [1 + /v cose D(/vcos)] (E. 13) N ( Zvcose ) J This last form is plotted in Fig. 6.9.

APPENDIX F THE MAGNITUDE OF THE PHASE RATE FILTER TRANSFER FUNCTION In finding the magnitude of the phase rate filter transfer function, the filter phase may be ignored. Without loss of generality, the impulse response of the filter may thus be centered about t = 0, and the mean correction terms in Eq. 6.4 may be ignored. Letting 0(T) be the 0 value at fixed f and time T, and including the time step A = 99.6 seconds explicitly, Eq. 6.4 becomes n 2 kA * 0(T+kA) -n 0(T) = (F.1) n 2 E (kA)2 -n where all summations are on k. Let N = - and ignore the scaling performed by the denominator. Then N (T) = E kA 0(T+k) (F.2) -N 182

183 Recognize by rewriting Eq. F.2 that this is simply the discrete convolution of 0 with an impulse response h(t) which looks like h(t) NA -NA N 0(T) = E -kA * N (T-kA) (F.3) -N The transfer function H(f) of the phase rate filter is the Fourier transform of h(t) N H(f) = -kA ~ ec, = 2Tf (F.4) -N To compute H(f), first compute G(f) N G(f) = e-jkA (F.5) -N G(f) is the sum of a geometric series. Using the wellknown formula for the sum of such a series,

184 jwNA -jw (N+1)A le -e G(f) = e -e (F.6). A -2 Factoring e from numerator and denominator yields -3 2 jw(N+ )A jw(N+ )A e e -e G(f) (F.7).wA A A (F.7) e-3 2 e 2 -e 2 Recognizing that ejx -e]x = 2j sin x sin('f(2N+1)A) G(f) (F.8) sin (rffA) Differentiating G with respect to f in Eq. F.5, ~dG(f)~ =f -j27 Tr k e- (F.9) -N or N -j dG(f) = -kA ewk = H(f) (F.10) 2~ df -N In terms of magnitudes dG|H(f) |||(F.11) II H(f) i 2~1 dGf~

185 Letting M = 2N + 1 be the number of e points used in the filtering operation and performing the differentiation of Eq. F.11, A Msin( TfA)cos(MTrfA.) - cos (rfA)sin (MTrfA) IH(f) o sin2 (TfA) (F.12) With no loss of applicability, M may be allowed to be even as this will affect the phase of H(f) but not the magnitude.

REFERENCES 1. Urick, R. J., Principles of Underwater Sound for Engineers, McGraw-Hill Book Company, New York, 1967. 2. Kailath, T., "Measurements of Time-Variant Comrmunication Channels," IRE Transactions on Information Theory. IT8, S229-S236, September, 1962. 3. Veenkant, R. L., Investigation of the Propagation Stability of -a Time Spread Underwater Acoustic Channel, Technical Report 226, Cooley Electronics Laboratory, The University of Michigan, Ann Arbor, May, 1974. 4. Birdsall, T. G., R. M. Heitmeyer, and K. Metzger, Modulation by Linear Maximal Shift Register Sequences: Amplitude, Biphase, and Complement-Phase -Modulation, Technical Report 216, Cooley Electronics Laboratory, The University of Michigan, Ann Arbor, December, 1971. 5. Van Trees, H. L., Detection, Estimation, and Modulation Theory, Part III, John Wiley and Sons, New York, 1971. 6. Heitmeyer, R. M., Underwater Sound Propagation in the Straits of Florida- The Preliminary Analysis of the MIMI Experiment of 1970, Technical Report 213, Cooley Electronics Laboratory, The University of Michigan, Ann Arbor, February, 1972. 7. Cederquist, G. N., A MIMI Propagation Study: Coherent Spectra of Wideband Underwater Acoustic Receptions in the Straits of Florida, 25 November 1970, Technical Report 215, Cooley Electronics Laboratory, The University of Michigan, Ann Arbor, May, 1973. 8. Unger, R. and R. Veenkant, Underwater Sound Propagation in the Straits of Florida: The MIMI Experiment of 3 and 4 February 1965, Technical Report 183, Cooley Electronics Laboratory, The University of Michigan, Ann Arbor, May, 1967. 9. Unger, R. and R. Veenkant, Underwater Sound Propagation in the Straits of Florida. The MIMI Continuous and Sampled Receptions of 11, 12, and 13 August 1966, Technical Report 186, Cooley Electronics Laboratory, The University of Michigan, Ann Arbor, June, 1967. 186

187 REFERENCES (Cont.) 10. Steinberg, J. C. and T. G. Birdsall, "'Underwater Sound Propagation in the Straits of Florida," Journal of the Acoustical Society of America, Vol. 39, No. 2, February, 1966, pp. 301-315. 11. Martin, J., Design of Man-Computer Dialogues, Prentice Hall, Inc., Englewood Cliffs, 1973. 12. Newman, W. M. and R. F. Sproull, Principles of Interactive Computer Graphics, McGraw-Hill Book Company, New York, 1973. 13. Williamson, H., "Algorithm 420: Hidden Line Plotting Program," Communications of the ACII, Vol. 15, No. 2, February, 1972, pp. 100-103. 14. Walters, R. F., "Contouring by I:lachine: A User's Guide," The American Association ofL Petroleum Geo." logists Bulletin, Vol. 53, No. 11_ November, 1969, pp. 2324-2340. 15. Akima., H., "A Method of Bivariate Interpolation and Fitting Based on Smooth Surface Local Procedures," Communications of the ACM, Vol. 17; No. 1, January, 1974, pp. 18-19 and 26-31. 16. Coulthard, W. J., UBC Pers-Perspective Views of Data on a Regular Grid, Computing Centre, University of British Columbia, Vancouver, Canada, October, 1969. 17. Hartwig, G. W., CONTUR. —A Fortran IV Subroutine for the Plotting of Contour Lines, Ballistic Research Laboratories Memorandum Report iNoe 2282, Aberdeen Proving Ground, Maryland, March, 1973 (NTIS Accession Number AD-760 437). 18. Hoare, C. A. R., "An Axiomatic Approach to Computer Programming," Communications of the ACM, Vol. 12, No. 10, October, 1969. 19. London, R. L., "Proving Programs Correct. Some Techniques and Examples," B.I.T., Vol. 10, 1970, pp. 168-172. 20. Dijkstra, E. S., "Structured Programming," in Structured Programming, Academic Press, 1972.

188 REFERENCES (Cont.) 21. Parnas, D. L., "On the Criteria To Be Used in Decomposing Systems into Mrodules," Communications of the ACM, Vol. 15, No. 12, December, 1972. 22. Parnas, D. L. and D. P. Siewiorek, Use of the Concept of Transparency in the Design of Hierarchically Structured Systems, Department of Computer Science, Carnegie-Mellon University, Pittsburgh, Pennsylvania, September, 1972. 23. Spier, M. J., "A System Theoretic Look at the Complexity of Access Control Mechanisms,' Proceedings of the International Workshop on Protection in Operating Systems, IRIA-LABORIA, Paris, France, August 13-14, 1974. 24. Wirth, N., Revised Report on the Programming Language PASCAL, Computer Science Department, Eidgentssische Technische Hochschule, Zurich, Switzerland, 1971. 25. Wirth, NEo, Systematic Programming, Prentice —Hall, Englewood Cliffs, New Jersey, 1973. 26. Neisser, U., "The Processes of Vision," Scientific American, Vol. 219, No. 3, September, 1968, pp. 204214.

DISTRIBUTION LIST Office of Naval Research (Code 412) 2 (Code 102-OS) 1 (Code 480) 1 Department of the Navy Arlington, Virginia 22217 Director 6 Naval Research Laboratory Technical Information Division Washington, D. C. 20375 Director Office of Naval Research Branch Office 1030 East Green Street Pasadena, California 91106 Office of Naval Research San Francisco Area Office 760 Market Street - Room 447 San Francisco, California 94102 Director Office of Naval Research Branch Office 495 Summer Street Boston, Massachusetts 02210 Office of Naval Research New York Area Office 207 West 24th Street New York, New York 10011 Director Office of Naval Research Branch Office 536 South Clark Street Chicago, Illinois 60605 Commanding Officer 8 Office of Naval Research Branch Office Box 39 FPO New York 09510 Commander 1 Naval Surface Weapons Center Acoustics Division White Oak Silver Spring, Maryland 20910 ATTN: Dr. Zaka Slawsky 189

190 DISTRIBUTION LIST (Cont.) Officer in Charge 1 Annapolis Laboratory Naval Ship Research and Development Center Annapolis, Maryland 21402 Commander Naval Sea Systems Command Code SEA 037 Washington, D. C. 20362 Commander Naval Sea Systems Command Washington, D. C. 20362 ATTN: Mr. Carey D. Smith (Code SEA 06H1) CDR Bruce Gilchrist (Code SEA 06H2) 1 Office in Charge 1 Pasadena Laboratory Naval Undersea Center 3202 East Foothill Boulevard Pasadena, California 91107 Commanding Of~ficer 1 Fleet Numerical Weather Central Monterey, California 93940 Defense Documentation Center 5 Cameron Station Alexandria, Virginia 22314 Chief of Naval Material 1 Department of the Navy Washington, D. C. 20360 ATTN: Mr. James Probus (Acting Director of Navy Laboratories) Commander Naval Electronic Systems Command Washington, D. C. 20360 ATTN: CAPT. J. Bajus (NAVELEX 03) 1 CDR. A. Miller (NAVELEX 320) 1 Chief of Naval Operations 1 Department of the Navy Pentagon, Room 5B718 Washington, D. C. 20350 ATTN: CAPT. Robert B. Brunsted

191 DISTRIBUTION LIST (Cont.) Commander 1 Naval Ship Research and Development Center Bethesda, Maryland 20034 ATTN: Mr. Craig Olson Chief of Naval Operations 1 Department of the Navy Pentagon, Room 4C559 Washington, D. C. 20350 ATTN: Capt. A. H. Gilmore Commander Naval Undersea Center San Diego, California 92132 ATTN: Dr. Dan Andrews 1 Mr. Henry Aurand 1 Chief Scientist 1 Navy Underwater Sound Reference Division P. O. Box 8337 Orlando, Florida 32806 Officer in Gharge 1 New London Laboratory Naval Underwater Systems Center New London, Connecticut 06320 Commander 1 Naval Air Development Center Warminster, Pennsylvania 18974 Commander 1 Naval Ship Research and Development Center Bethesda, Maryland 20084 Superintendent 1 Naval Postgraduate School Monterey, California 93940 Commanding Officer 1 Naval Coastal Systems Laboratory Panama City, Florida 32401 Commanding Officer 1 Naval Underwater Systems Center Newport, Rhode Island 02840 Superintendent 1 Naval Academy Annapolis, Maryland 21402

192 DISTRIBUTION LIST (Cont.) Commanding Officer Naval Intelligence Support Center 4301 Suitland Road Washington, D. C. 20390 ATTN: Dr. Johann Martinek Mr. E. Bissett Commander Naval Sea Systems Command Code SEA 03E Washington, D. C. 20362 Dr. Melvin J. Jacobson Rensselaer Polytechnic Institute Troy, New York 12181 Dr. Charles Stutt General Electric Company P. O. Box 1088 Schenectady, New York 12301 Dr. Alan Winder MSB Systems, Inc. 110-16 72nd Avenue Forest Hills, New York 11375 Dr. T. G. Birdsall Cooley Electronics Laboratory University of Michigan Ann Arbor, Michigan 48105 Dr. Harry DeFerrari University of Miami Rosenstiel School of Marine and Atmospheric Sciences Miami, Florida 33149 Mr. Robert Cunningham Bendix Electronics Center 15825 Roxford Street Sylmar, California 91342 Dr. Stephen Wolff John Hopkins University Baltimore, Maryland 21218 Dr. M. A. Bas in S. D. P., Inc. 15250 Ventura Boulevard, Suite 518 Sherman Oaks, California 91403

193 DISTRIBUTION LIST (Cont.) Commanding Officer New London Laboratory Naval Underwater Systems Center New London, Connecticut 06320 ATTN: Dr. Albert Nuttall Dr. Walter Duing University of Miami Rosenstiel School of Marine and Atmospheric Sciences Miami, Florida 33149 Commanding Officer New London Laboratory Naval Underwater Systems Center New London, Connecticut 06320 ATTN: Dr. H. W. Marsh Dr. D. M. Viccione Dr. David Middleton 127 East 91st Street New York, New York 10028 Dr. Donald W. Tufts University of Rhode Island Kinston, Rhode Island 02881 Dr. Loren W. No te Department of Electrical Engineering FT-10 University of Washington Seattle, Washington 98195 Mr. S. W. Autrey Hughes Aircraft Company P. O. Box 3310 Fullerton, California 92634 Dr. Thomas W. Ellis Texas Instruments, Inc. 13500 North Central Expressway Dallas, Texas 75231 Mr. Rohert Swarts Applied Physics Laboratory University of Washington 1013 Northeast Fortieth Street Seattle, Washington 98195

194 DISTRIBUTION LIST (Cont.) Institute for Acoustical Research 3 Miami Division of the Palisades Geophysical Institute 615 S. W. 2nd Avenue Miami, Florida 33130 ATTN: Mr. M. Kronengold Dr. J. Clark Dr. C. Kimball Mr. Carl Hartdegen Palisades Geophysical Institute Sofar Station FPO New York 09560 Mr. Charles Loda Institute for Defense Analyses 400 Army-Navy Drive Arlington, Virginia 22202 Mr. Beaumont Buck Polar Research Laboratory 123 Santa Barbara Avenue Santa Barbara, California 93101 Dr. M. Weinstein Underwater Systems, Inc. 8121 Georgia Avenue Silver Spring, Maryland 20910 Dr. Thomas G. Kincaid General Electric Company P. O. Box 1088 Schenectady, New York 12301 Applied Research Laboratories 4 The University of Texas at Austin P. O. Box 4029 Austin, Texas 78712 ATTN: Dr. Lloyd Hampton Dr. Charles Wood Dr. T. D. Plemons Woods Hole Oceanographic Institute Woods Hole, Massachusetts 02543 ATTN: Dr. Paul McElroy Mr. R. Porter Mr. R. Spindel

195 DISTRIBUTION LIST (Cont.) Dr. John Bouyoucos Hydroacoustics, Inc. 321 Northland Avenue P. O. Box 3818 Rochester, New York 14610 Systems Control, Inc. 1 260 Sheridan Avenue Palo Alto, California 94306 ATTN: Mr. L. Seidman Atlantic Oceanographic and Meteorological Laboratorie s 15 Rickenbacker Causeway Miami, Florida 33149 ATTN: Dr. John Proni Dr. C. N. K. Mooers University of Miami Rosenstiel School of Marine and Atmospheric Sciences 10 Rickenbacker Causeway Miami, Florida 33149 Westinghouse Electric Corporation Advanced Development Programs Marketing Department -MS 227 P. O. Box 746 Baltimore, Maryland 21203 ATTN: F. J. Frissyn Oak Ridge National Laboratory 1 Union Carbide Corporation Nuclear Division P. O. Box X Oak Ridge, Tennessee 37830 Cooley Electronics Laboratory 100 University of Michigan Ann Arbor, Michigan 48105

UNIVERSITY OF MICHIGAN 3 9015 02947 4916 3 9015 02947 4916