DEPARTMENT OF ELECTRICAL ENGINEERING 

SCHOOL OF ENGINEERING 

OLD DOMINION UNIVERSITY 

NORFOLK, VIRGINIA 

- 

0 

rr) 

M 

r* 

vl 

;d -0 

I pi 23 

iB 

DESIGN OF INFRASOUND-DETECTION SYSTEM 

us 

a a Ul 

VIA ADAPTIVE LMSTDE ALGORITHM 2i 3 - 

r7 

c3 

\ 

9 

'3 

By - .J 

Q 

Camille S. Khalaf -4 

2- *b, 

u 15 s r-Q 

and R Q -4 "J n( -4 3 

4 543 

John W. Stoughton, Principal Investigator Q 4 ;U a3 

0 a 

Final Report ~Z.LO Orzl~~ For the period ending December 31, 1984 i5-4 -t~ 

=v) 21 aa*-#r-i 

'-'LCaa 

rn aFi7 = -4 

QOrySV 

n wo 

Prepared for the -'J ~-1 .-- ma 

National Aeronautics and Space Administration ~.am k-4 M . z Langl ey Research Center m rq i-i u \. 

Hampton, Virginia 23665 p1 QQP;Q~- 1399 

;&adcqt.J 

3 0 x 

I mwp Under ~MQG

Cr) D3R.4 Research Contract NAS1-17099 4 ik C,Q qg a 

Task Authorization No. 33 -~+a~t- ===.=rlycr, 

Dr. A1 lan J. Zuckerwar, Technical Monitor 

IRD-Acoustics & Vibration Instrumentation Section 

November 1984 
OEPAKTMENT OF ELECTRICAL EhGINEERING 

SCHOOL OF ENGINEEKIWG 

ULD DOMINION UWIVERSITY 

NOKFOL K, VIRGINIA 

DESIGN OF INFRASOUND-DETECTION SYSTEM 

VIA ADAPTIVE LMSTDE ALGORITHM 

~amille' S. Khalaf 

and 

John A. Stoughton, Principal Investigator 

Final Report 

For the period ending December 31, 1984 

Prepared for tne 

National Aeronautics and Space Administration 

Langley Research Center 

Hampton, Virginia 23665 

Under 

Research Contract NAS1-17099 

Task Authorization No. 33 

Dr. A1 lan J. Zuckerwar, Technical Monitor 

IRD-Acoustics & Vibration Instrumentation Section 

Submitted by the 

01 d Dominion University Research Foundation 

P. U. Box 6369 

Norfol~, Vi rginia 23508 

November 1984 
ABSTRACT 

DESIGN OF INFRASOUND-DETECTION SYSTEIq 

VIA AUAPTIVE LMSTDE ALGORITHM 

Camille S. Khalaf 

01 d Dominion University 

Di rector: John W. Stoug hton 

A proposed solution to an aviation safety problem is based on pass i ve detection of turbulent weather phenomena through their i nfrasoni c 

emission. This thesis describes a system design that is adequate for 

detection and bearing evaluation of infrasounds. An array of four sensors, with the appropriate hardware, is used for the detection part. 

Bearing evaluation is based on estimates of time delays between sensor 

outputs. The general i zed cross correl ation (GCC) , as the conventional 

t ime-del ay estimation (TDE) method, i s first reviewed. An adaptive TOE 

approach, using tne least mean square (LiqS) algorithm, is then discussed. A comparison between the two techniques is made and the advantages 

of the adaptive approacn are listed. The behavior of the GCC, as a Koth 

processor, is examined for the anticipated signals. It is shown that 

the Roth processor has the desired effect of sharpening the peak of the 

correlation function. It is also shown that the LMSTDE technique is an 

equivalent implementation of the Roth processor in the time domain. A 

LMSTDE lead-lag model, with a variable stability coefficient and a convergence criterion, is designed. This model is employed in an automatic 

scheme developed for the sensor array. The software and hardware system 

parameters are derived and determined. The effectiveness of the system 

is illustrated through simulation and field testing. 
ACKNOWLEDGEMENTS 

I wish to thank my thesis advisor, Dr. John W. Stoughton, for his 

invaluable assistance and insights, and for his support and understanding when I needed them most. I wish to acknowl edge my gratitude to Dr. 

A1 1 an J. Zuckerwar, Acoustical and Vi brational In strunentation Section 

(AVIS), NASA/Langley Research Center, whose constant effort and assistance made this research possible. My thanks also go to Mr. Harlan K. 

Holmes, Head-AVIS, Mr. Ralph Lewis, Aerospace Technician-AVIS, and Mr. 

David Katzoff, Governor's School, for their time and effort. I acknow1 edge the financial support by the National Aeronautics and Space Administration under contract NAS1-17099, Tasks 10 and 33. 
TABLE OF CONTENTS 

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

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

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

UEFINITION OF SYMBOLS .......................................... 

1 . INTRODUCTION ........................................ 

1.1 Motivation ..................................... 1.2 Turbulence Detection ........................... 1.3 Research Objectives ............................ 

2 . THEORY ................................................ 

2.1 Introduction ................................... 2.2 General i zed Cross Correl ation Approach ......... 2.3 Roth Processor ................................. 2.4 Adaptive Filters ............................... 2.5 Adaptive LMSTDE Approach ....................... 2.6 Summary ........................................ 

3 . SYSTEM DESIGN ....................................... 

3.1 Introduction ................................... 3.2 Design Assumptions ............................. 3.3 Filter Design .................................. 3.4 Specifications and Hardware .................... 3.5 Summary ........................................ 

4 . SIMULATION AND RESULTS .............................. 

4.1 Introduction ................................... 4.2 System behavior ................................ 

4.2.1 Nominal Response ........................ 4.2.2 Filter Parameters'....................... 

4.2.3 Signal Parameters ....................... 4.2.4 Shuttle Launch Noise .................... 

4.3 System Operati on and Results ................... 4.4 Man-Made Signals ............................... 

. 

PAGE 

i i 

iii 
PAGE 

4.5 Weather-Re1 ated Data ........................... 102 

4.6 Summary ........................................ 115 

5.1 Remarks ....................................*... llb 

5.2 Future dork .................................... 118 

REFERENCES .................................................... 

APPENDIX A .................................................... 
- TABLE 

LIST OF TABLES 

- PAGE 

2.1 Correl ation weighting functions.. . . . . . . . . . . . . . . . . . . . . 14 

3.1 Parameter specifications ............................. 50 
LIST OF FIGURES 

FIGURE 

1.1 

PAGE 

Early warning system .................................. 

klay-to-angle conversion............................. 

Generalized cross correlation. ........................ , 

Adaptive 1 inear combinatori a1 system.. ................ 

Modeling TDE using adaptive filter .................... 

TDE lead-lag model. ................................... 

LMSTDE a1 gori t hm' s flowchart. ......................... 

Variance computation .................................. 

Delay estimate computation ............................ 

System block diagram .................................. 

Type of random data sequences used for s imul at i on. .... 

Nominal filter response for D = +4 samples (C=10, 

N=17, L=128, SNR=l). .................................. 

Nominal filter res onse for D = -4 samples (C=10, N=17, L=128, sNR=lP.. ................................. 

N=7. Delay cannot be estimated with N<2xDmax+l 

(D=4, C=10, L=128, SNR=l). ............................ 

N=33. No improvement over N=17. (D=4, C=10, 

L=128, SNR=l) ......................................... 

Instabi 1 i ty for C=0.5. (D=4, N=17, L=128, SNR=l). .... 

C=200. No improvement over C=10. (D=4, N=17, 

L=128, SNR=l) ......................................... 

L=32. hbiguity in filter output due to short 

record length. (D=4, C=10, N=17, SNR=l). ............. 
FIGURE 

4.9 

- PAGE 

70 

L=512. Improved response over L=128 but not 

necessary. (D=4, C=10, N=17, SNR=1) .................. 

SNR=1/2. Sl ow convergence. (D=4, C=10, N=17, 

L=128) ................................................ 

SNR=2. Faster convergence. (D=4, C=10, N=17, 

L=128) ................................................ 

7r 

WI=~. Main lobe width of 8 rdso~utions. (D=4, 

C=lO, N=17, L=128, SNR=l). ............................ 'li 

wl =x. Main lobe width of 16 resolutions. (D=4, 

c=lO, N=17, ~~128, SNR-1). .............. ..: ........... 

Filter output for time varying delay (1 point 

per 64 points) ........................................ 

Filter out ut for time varying delay (1 point per 128 points ! ........................................... 

Array configuration used in sampl i ng STS-8 1 aunch 

signal ................................................ 

Time histories of STS-8 launch signal (in equivalent 

arbitrary units) on one channel.. ..................... 

Auto power spectrums of STS-8 1 aunch signal 

(in arbitrary units) on one channel ................... 

Time hi stories (in equivalent arbitrary units) of 

STS-8 launch signal on three microphones. ............. 

Filter response for launch signal, C2 x C3.. .......... 

Filter response for 1 aunch signal, C3 x C4.. .......... 

Filter response for launch signal, C4 x C2.. .......... 

True direction of shuttle launch site.. ............... 

Software f1 owchart.. .................................. 

Location of the runway in Langley Air Force Base 

relative to the microphone array............. ......... 
FIGURE 

Time histories (in equivalent arbitrary units) of 

i nfrasounds generated by jet powered aircraft. ........ 

Auto-power spectrums (in arbitrary units) of 

infrasounds generated by jet powered aircraft.. ....... 

Time histories (in equivalent arbitrary units) 

of infrasounds generated by jet powered aircraft 

sampled on (a) channel 1, (b) channel 2, and 

(c) channel 3.. ....................................... 

Filter response for jet data, C1 x C2 ................. 

Filter response for jet data, C2 x C3 ................. 

Filter response for jet data, C3 x Cl.. ............... 

Time hi stories (in equivalent arbitrary units) of hurricane signal sampled on one channel at 

different times ....................................... 

Auto-power spectrums (in arbitrary unit s) of 

hurricane signal sampled on one channel at 

different times.... ................................... 

Auto-power spectrums on a log scale (log amplitude 

vs. log frequency, in equivalent arbitrary units) 

of hurricane signal sampled one channel at different 

times ................................................. 

Time hi stories (in equivalent arbitrary units) 

of hurricane signal simul taneousl y sampl ed at 

(a) channel 1, (b) channel 2, (c) channel 3 ........... 

Auto-power spectrums (in arbitrary units) of 

hurricane si nal simul taneousl sampled at 

(a) channel 7, (b) channel 2, fc) channel 3.. ......... 

Auto-power spectrums on a 1 og scal e (1 og-amp1 itude 

vs. log-frequency, in equivalent aribtrary units) 

of hurricane signal simultaneously sampled at 

(a) channel 1, (b) channel 2, (c) channel 3. .......... 

Filter response for hurricane data, C1 x C2 ........... 

Filter response for hurricane data, C2 x C3 ........... 

Filter response for hurricane data, C3 x C1 ........... 

PAGE 

96 

9 7 

viii 
LIST OF FIGURES - Concluded 

Figure Page 

4.41 Hurricane site re1 ative to the microphone array 

inNASALangleyResearchCenter ....................... 114 

A.l Elevationeffectondelays ............................ 123 

A.2 Re1 ative angles produced by a hori zontal wave. ........ 125 
DEFINITION OF SYMBOLS 

continuous time 

discrete time 

time 1 ag 

time de1 ay 

frequency in Hz 

sampl i ng rate 

sampl ing period 

r ad i an frequency 

di screte del ay 

bearing re1 ative to microphones axis 

bearing relative to North 

elevation angle 

bearing resol ution 

distance between two sensors 

sound ve1 oci ty 

stability constant 

scale of stabil ity constant 

filter length 

record length 

file length 

variance window 

variance threshold 
signal 

noi se 

weighting 

delta function 

input signal 

random sequence 

fi 1 ter error function 

filter desired response 

f i 1 ter output 

filter input vector 

correl ation function 

power spectrum 

cov ari ance matrix 

coherence function 

Fourier transform 

expected v a1 ue 

gradient 

convol uti on 
CHAPTER 1 

I NTROD UCT I0 N 

1.1 Motivation 

Current research programs have lead to significant advances in 

ground-based and airborne equipment for providing i nformati on re1 ative 

to severe turbulent weather. However, most of these programs are still 

in their experimental stages and very few are operational . There continues to De a serious aviation safety probl€?m associated witn aircraft 

operations in the vicinity of severe storms. 

In 1981, general aviation aircraft numbered more tnan 200,000 and 

flew more than 40 mill ion hours [I]. General aviation operations 

accounted for 662 fatal accidents from all causes, with 1,265 fatalities 

(FAA, 1981). Informal accident cause/f actor statistics from the National Transportation Safety Board for 1981 indicate that weather caused, or 

was a factor, in 40 percent (289 cases) of the U. S. general aviation 

accidents. Earlier statistics indicated that turbulence is the largest 

single cause of weather-related air carrier accidents in the U. S. From 

1962 to 1974, turbulence was either a cause of or a contributing factor 

in 189 of 450 weather-related cases [2]. 

The two types of turbulence usually encountered are clear-air tur- 

~ul ence (CAT) and thunderstorm-re1 ated turbulence. CAT, a problem for 

all aircraft, cannot be seen because it usually has no cloud signature. 

It may develop in a standing wave caused by air noving over mountainous 
terrain, and is frequently associated with shear-induced KelvinHelmholtz atmospheric waves occurring in a stati stical ly stable atmosphere [3]. Accidents caused by CAT are not as serious as the ones 

related to turoulence associated with thunderstorms. A CAT accident may 

result in discomfort, injuries, aircraft damage, and/or unscnedul ed 

1 andi ng. Thunderstorms and other convective clouds are critical 1 y 

important sources of lowral titude turbulence and wind variabi 1 ity. Many 

produce strong downdrafts that transport air downward, which then 

spreads out rapidly over the ground. This mechanism, if encountered 

during take off or landing of aircraft, may result in serious if not 

fatal accidents. 

One of the pressing aviation safety problems is that of providing 

the pilot with information needed to avoid turbulence hazards which 

exceed the design capabilities of the airplane. A proposed solution to 

this proolem through passive detection of turbulence i s presented in the 

next section. 

1.2 Turbulence Detection 

The primary technique for detecting turbulence and storm cells is 

Doppler radar. This technique requires tne presence of reflective 

particles such as precipitation or dust for meter01 ogi cal appl i cations. 

Doppler radar proved ineffective in cases of CAT or developing storm 

cells due to the absence of reflective particles in these phenoiena. 

The proposed technique consists of passive detection of 1 arge-scale 

patches of turbulence i n the Earth's atmosphere through infrasonic 

emissions. The infrasonic technique offers the advantages of being (1) 
passive, (2) inexpensive, and (3) inherently more sensitive to atmospheric di sturbances than electromagnetic or optical techniques. The 

assuned utilization of the infrasonic technology is illustrated in Figure 1.1. Infrasonic emissions from the patch of turbulence are detected 

at two stations, each containing an array of four sensors. By means of 

real-time signal processing, the direction of the turbulent source, 

based on time delays between sensors, is determined at each station. 

The turbulence is then located by triangulation between the two stations. This information is relayed from a control center to a flight 

services advisory, and from there to the pilots of approaching aircraft. 

A series of such stations would provide early warning along the length 

of domestic air traffic routes. 

It is informative at this point to mention that the first step 

toward realization of the infrasonic technique was the development of a 

unified acquisition system for acoustic data [4] by Dr. A1 lan J. 

Zuckerwar and Mr. Harl an K. Holmes at NASA Langley Research Center. 

This system served an important role in understanding the detection 

problem during the early stage of this research. 

The del ay-to-angle conversion needed in this passive approach i s 

illustrated in Figure 1.2. If rd is the time delay between the 

arrival of the infrasonic signal at two spatially separated sensors, 

then the direction of the turbulent source is given by 
Figure 1.1. Early warning system. 
Figure 1.2. Del ay-to-angle conversion. 
where V is the sound velocity (1080 ft/sec at O°C), and 2 is the 

distance between the two sensors. For a given delay, the solution of 

equation (1.1) can be interpreted in two ways as shown in the figure. 

The signal associated with the false direction is usually referred to as 

the "ghost" signal. The presence of the ghost-signal requires a minimum 

of three sensors in the array to uniquely determine the source direction. 

1.3 Research Objectives 

The infrasonic technique, proposed to provide pilots with an early 

warning of turbulence, was described in the previous section. The three 

aspects that are essential to its success are the following: 

1. emission and propagation of infrasounds from turbulent weather 

phenomena; 

2. detection of infrasounds and bearing evaluation of their 

sources at every station; 

3. locating the sources by triangulation between stations via a 

communication network that is a1 so responsible for transferring 

the information to pilots. 

The first aspect is out of the designer's control and is based on the 

theoretical predictions of infrasonic emission and propagation. The 

third aspect constitutes the commercial phase of the program once the 

performance of the individual stations proves satisfactory. 

The goal of this research is the overall system design that is 

adequate for imp1 ementi ng tne second aspect of the infrasonic technique, 

In this context, the system should function as an infrasonic detector, 
time-del ay estimator, and del ay-to-angl e converter. The detection part 

i nvol ves the sensor array and correspondi ng hardware b1 ocks wh i 1 e the 

del ay-to-angle conversion is accompl ished by a simple software routine. 

The main design issue is that of an appropriate signal processing tecnnique for time-delay estimation. The GCC method and the LMS parametric 

technique will both be reviewed. More specifically, the GCC method 

using a Roth weighting function, and the adaptive implementation of the 

LMS technique, will be examined. Design decisions, regarding the TDE 

technique as well as the hardware system, wi 11 then be made according to 

the anticipated signal characteristics. 

In Chapter 2, the theory behind TDE techniques is reviewed. Based 

on this theory and the anticipated infrasounds, the system design is 

presented in Chapter 3. Chapter 4 evaluates the system through simulation of the basic TDE algorithm response and discussion of the system 

performance in the actual field. Conclusions and future research are 

presented in Chapter 5. 
CHAPTER 2 

THEORY 

2.1 Introduction 

The problem identified in chapter one is a fundamental passive 

sonar signal processing problem in which delays between the times of 

arrival of the pertinent acoustic waves at four sensors are to De est imated. This chapter will discuss the theory behind TDE techniques. 

Section two reviews the conventional generali zed cross correl ation 

approach (GCC) while section three reviews the Rotn processor specif ically. Uf particular interest is a parametric approach through Widrow's 

adaptive filter. The filter structure is reviewed in section four while 

its application to TOE is presented and thoroughly investigated in the 

fifth section. The chapter is concluded in section six by presenting 

the advantages of the adaptive least mean squared time de1 ay estimation 

(LMSTDE) method. 

2.2 General i zed Cross Correlation Approach 

The GCC approach to time delay estimation has been discussed by 

many i nvest igators . We1 1 known references i nc1 ude papers written by 

Knapp and Carter [5], and Hassab and Boucher [6]. This section only 

reviews this approach through Figure 2.1. xl (t) and x2(t) are sampled at two spatially separated sensors and then fed to a basic cross 

correl ator. The ~asic cross correlator, as discussed by Papoul is [7], 
Inverse 

Fourier 

Transform 

I 

I I 

I 0 x2(t) 0 Xl(t) , 

Sensors 

I *denotes 

I complex 

I conjugate. 

Linear 

Fi l ter 

------ --1 

---- 1 

L--- 

r----'- 

I Peak 

---- ----- --- - < 

___----- 

I Detector I 

I V 

Time delay estimate 

v I 

Figure 2.1. Generalized cross correlation. 

I 

I 

I 

Fourier 

Transform 

I 

I 

I 

Fourier 

Transform 

I I 

I 

I 

I 

I 

I * I Basic i 

I I cross 

I correlator 
computes the cross correlation function, R (T) , between xl (t) and 

x1 x2 

X* (t) by means of the inverse Fourier transform of their cross po$r 

spectrum, G (f) . Then the delay estimate is simply the abscissa 

x1 x2 

value at which the cross correlation function peaks. In order to improve the accuracy of the delay estimate, a linear filter, w(t), is 

convolved with tne output of the cross correlator. A peak detector 

routine is then used to determine the abscissa value of the peak in the 

filtered correlation function. In practice, a frequency weighting d(f) 

= F [w(t)], where F[*j denotes the Fourier transform of [*J, which is 

equivalent to w(t) , is appl ied to the cross power spectrum prior to 

taking the inverse Fourier transform. This frequency weighting replaces 

the linear filter so that a11 computations, except for the peak detector 

routine, are done in the frequency domain. A discussion of this 

weighting, after mathematically modeling the system, is presented next. 

A signal emanating from a remote source and monitored in the presence of noise at two spati ally separated sensors can be mathematically 

modeled as 

xl (t) = Sl (t) + ml(t) 

where Sl (t), ml (t) and rn2(t) are real, zero-mean, jointly stationary 

random processes, rd denotes time delay and a an attenuation factor, 

Signal Sl (t) is assumed to be uncorrel ated with ml (t) and m2 (t) . 
The cross correlation between xl (t) and x2(t) is related to the 

cross power spectral densityfunction, G (f), bythe Fouriertransx1 x2 

form re1 at ionshi p 

w 

R (T) = 1 G (f) e jaf~ df

X1 X2 a x1 x2 

In practice, only an estimate G; (f) of (2.2) can be obtained from 

12 

finite observations of xl(t) and x2(t). Consequently, the output of 

the basic cross correlator, when no weighting is used, is: 

It is informative at this point to examine the shape of the cross 

correlation function R (r ), so that the necessity of using a x1 x2 

weiyhtingfunctioncanbejustified. Anexpressionfor R (T) can 

x1 x2 

be obtained from (2.1), using the expectation operator €[a], as 

The Fourier transform of (2.4) gives the cross power spectrum 

-jmf-c 

G (f) = aG (f) e + G (f) . x1 x2 Sl S1 mi m2 
If ml (t) and % (t) are uncorrel ated, 

Gml m2 

(f) = 0, and since i 

multiplication in one domain is a convolution in the transformed domain, 

it follows from (2.5) that 

-j&f~ 

wnere B denotes convolution and 6 (T) = F- ie me 

interpretation of (2.6) is that the delta function has been spread or 

"smeared" by the Fourier transform of the signal spectrum. If Sl (t) 

is a white noise source, then its autocorrelation function is a delta 

function and no spreading takes place. However, for most practical 

applications this is not the case and spreading acts to broaden the peak 

of the cross correlation. 

To minimize the spreading effect, many weighting functions have 

been proposed in the literature (see Table 2.1) to operate on the cross 

power spectrum given in (2.5). With a general weighting, dg(f), the 

estimate of the genera1 i zed cross correl ati on becomes 

w 

I@ (T) = J ~g(f) G (f) e jhf~ df 

x1 x2 -OD x1 x2 

Wg(f) should be chosen to ensure a 1 arge sharp peak in I?g (r) rather 

x1 x2 

than a broad one in order to obtain a good time-delay resolution. 

2.3 Rotn Processor 

The weighting functions found in the literature are listed in Table 
2.1, where the notation y(f) has been used for the coherence function 

which is defined as 

The selection of Wg(f) ,to optimize certain performance criteria has 

been studied by several investigators (see, for example [8]). The purpose of considering the Roth processor is that it is equivalent to the 

time domain weighting involved in the adaptive LMSTDE approach. This 

equivalence will be shown in section five and will be helpful in gaining 

insight to the adaptive method. 

The weighting proposed by Roth [9] is 

Substituting for Wg(f) in (2.7) yields 

Equation (2.8) estimates the impulse response of the optimun linear 

(Wiener-Hopf) filter [lo] 
Tab1 e 2.1. Weighting Functions 

Processor Name Weight Reference 

Cross Correl ati on 

Ro th Processor 

Eckart Fi1 ter 

Smoothed coherence transform 

Ph ase transform 

3Maxirnum 1 i kel i hood 
which "best" approximates the mapping of xl (t) to x2(t) . If ml (t) # 

0, then Gx , (f) = GS1sl(f) + Gqml (f), and equation (2.6) becomes for 11 

the Roth processor 

From (2.9) we can conclude that only when G (f) is negligiole or 

ml ml 

when it equals any constant times s1 (f), the spreading does not 

occur and R (T) becomes a delta function. However, as can be seen 

12 

from the integral part in (2.9), the Roth processor has the desira~le 

effect of suppressing those frequency regions where G, (f) is 1 arge 1 "'1 

and Gx (f) is more likely to be in error. 

12 

The general i zed cross correl ation approach suffers regardless of 

which weighting function is used, from two basic facts. First, it 

relies on a sequence of fast Fourier transform (FFT) computations that 

tends to be time consuming. Second, as seen in Table 2.1, it requires a 

priori knowledge of the signal and noise statistics to implement the 

specific weighting. In passive detection problems, this information is 

unknown and if it were to be estimated, it would increase the complexity 

of the process and the time involved. The next section presents the 

adaptive LMSTOE algorithm, which is equivalent to the Roth processor, 

and is able to overcome both difficulties found in the generalized cross 

correl ati on method. 
2.4 Adaptive Filters 

This section will review the basic structure of adaptive filters 

discussed by Widrow [ll], while section 2.5 presents the appl ication of 

the filter in TDE problems. 

A signal filtering approach using an adaptive filter is in some 

sense a self-designing (real 1y self-optimi zing) process. The adaptive 

f i 1 ter described here bases its own "design" (its internal adjustment 

settings) upon estimated or measured statistical characteristics of the 

input and output signals. The statistics are not measured expl icitly; 

rather, the filter design is accomplished in a single process by a recursive LMS algorithm that automati cal ly updates the system coefficients 

with the arrival of each set of data samples. Figure 2.2 illustrates 

schematically the adaptive filter used in this case as a 1 inear combinatorial system. The filter consists of a set of variable weights (filter 

coefficients) whose input are the sampled input signals, a summer to add 

the weighted signals, and an algorithm to adjust the weights automatically. The impul se response of such a discrete system is completely 

control led by the weight settings. The adaptation process automatical 1 y 

seeks an optimal filter impulse response by adjusting the weights using 

gradient techniques to minimize the mean-square-error function. 

The analysis of the adaptive filter can be developed by assuming 

that the input signal s are statistical ly stationary random processes. 

Let the nth set of input signals be a vector X(n) of length N, 

xT(n) = xp(n) *** 
d m aJ- - 

-rb s 

as 

E'n x 

rb 'r 

V) m 
where xT denotes the transpose of X. Let the set of weights, at the 

nth time, be designated by the vector wT(n) = [wl(n) wdn) . . . wN(n)]. 

N 

The nth output signal is y(n) = Z wi(n) xi (n). This can be written 

i =l 

in matrix form as 

T y(n) = W (n) x(n) = xT(n) W(n). (2.10) 

Denoting the desired response by d(n), the error at tne nth time is 

T e(n) = d(n) - y(n) = d(n) - W (n) X(n). (2.11) 

The square of tne error is 

c2 (n) = & (n) - 2d(n) xT(n) W(n) + wT(n) X(n) xT(n) W(n) . (2.12) 

The expected value of c2 (n) is 

T T E [c2 (n)] = $ (n) - 2 Vdx W(n) + W (n) Vxx W(n) . (2.13) 

where (T) denotes the mean value of (a ), or E[*]. The vector Vdx is 

the cross covariance between d(n) and x(n) and is defined as 
and VxX denotes the auto-covariance matrix of X(n) 

It may be observed from (2.13) that for stationary input signals, the 

mean-square error is a second-order function of the weights. Thus, the 

mean-square-error function may be viewed, as suggested by Widrow [ll] , 

as a "performance surface" for the adaptive process that has a unique 

stationary point (minimum) which can be sought using gradient techni ques. 

The gradient at any point on the performance surface can be obtained by differenti ating the mean-square-error function (2.13) with respect 

to the weights. The gradient is 

The "optimal " weight vector, WLJqS, that yields the least-mean-square 

error, is found where the gradient is zero. Accordingly, 
where [* 1-l denotes the inverse of [* 1. Equation (2.17) is the Wienerkpf equation in matrix form [12]. Solving (2.17) for the optimum 

weight vector (the minimum point on the performance surface) can present 

serious computational problems. However, it will be shown that the 

adaptation process tries to find an exact or an approximate solution to 

the Wiener-Hopf equation by using the LMS algorithm with less computati onal complexity. , 

When using the LMS algorithm, changes in the weight vector are made 

along the direction of the estimated gradient vector. Accordingly, 

W(n+l) = W(n) + Ks v [s (n)] (2.18) 

where 

W(n) & -weight vector before adaptation 

W(n+l) 4 - weight vector after adaptation 

K 4 scalar constant controlling the rate of convergence and 

S - 

stability (Ks < 0) 

- 

V' [g2 (n) estimate of gradient of E[E~] = ez with respect to W, 

with W = W(n) 

One method for obtaini ng the estimated gradient of the mean-square-error 

function is to take the gradient of a single time sample of the squared 

error. That is, 
From equation (2.11), 

T v[e(n)] = V [d(n) - W (n) X(n)] = - X(n) 

Thus, 

~[;2 (n) ] = - 2 E (n) X(n) (2.20) 

It can be shown that the gradient estimate of (2.20) is unbiased, so 

that 

E { v [;2(n) j = v [F2(n)] 

Substituting (2.20) in (2.18) yields 

and the next weight vector is obtained by adding to the present weight 

vector the input vector scaled by the value of the error. This is the 

LMS algorithm. 

Next, the expected value E[W(n)] of the weight vector after a 

large number of iterations will be shown to converge to the Wiener solution given by (2.17). For this purpose, assume that the time between 

successive iterations of the algori thrn i s sufficiently long so that tne 

input vectors X(n) and X(n+l) are uncorrelated. Taking the expected 
value of both sides of (2.21) we obtain a difference equation in the 

expected value of the weight vector: 

T E[W(n+l)] = E[W(n)] - 2 Ks E { X(n) [d(n) - X (n) ~t(n)]} 

where I is an identity matrix. With an initial weight vector 'd(U), 

n+l iterations of (2.22) yield 

Equation (2.23) can be put in diagonal form by using the normal-form 

expansion of the matrix V. That is, 

where the diagonal matrix of eigenval ues is A, and the square matrix 

of eigenvectors is the modal matrix Q. Equation (2.23) may now be 

expressed as 
Consider the diagonal matrix [I+2 KsA]. As long as its diagonal terms 

are all of a magnitude less than unity, 

rim [I+ 2 K~A]"' + 0 

n+- 

and the first term of (2.24) vanishes as the number of iterations .increases. The sunmation factor in the second term of (2.24) becomes 

where the formula for the sum of a geometric series has been used. That 

is, 

Thus, in the 1 imit, equation (2.24) becomes 

which is the Wiener-Hopf solution in (2.17). 

Convergence of E[W(n+l)] to (2.17) is obtained if and only if the 

diagonal terms of [I+2 KSA] all have magnitudes less than unity, and 

since all eigenvalues in A are positive (the auto-covariance matrix, 
"XX is positive definite), the bounds on KS are given by 

where Am is the maximum eigenvalue of Vxx. The convergence condition 

can be related to the total input power as follows: 

N 

Am ( trace [Vxx] = = Total input power. 

i =l 

Theref ore ; 

For a slow, precise adaptation KS is usually chosen such that 

It is believed that the assumption of independent successive input 

vectors used for the convergence proof is overly restrictive. Griffiths 

[13J has shown that adaptation using highly correl ated successive sampl es converges to the Wiener solution, but 1 eads to higher steady-state 

mean-square error. Thus, we can conclude that for short-term stationary 

signals, with the feedback estimate, Ksy being bounded by (2.25), the 
weight vector, after a large nunber of iterations, is expected to converge to the Weiner solution that best maps the input x(n) to the 

desired response d(n) . 

2.5 Adaptive LMSTDE Approach 

The LMS adaptive filter described in the previous section has 

been widely applied in situations where the statistics of the inputs are 

either unknown or parti a1 ly unknown. Some appl i cat ions incl ude noise 

cancel 1 ing [14], 1 ine enhancing [14], prediction [Ill, spectrum analysis 

[13], and adaptive array processing [15]. A recent application of the 

filter is time delay estimation by F. Reed, P. Feintuch and N. Bershad 

[16]. This application is demonstrated in Figure 2.3, where the adaptive filter has a slightly different structure. The filter to be considered here has only two inputs: a primary input ~(n) and a secondary input (desired response) d(n) . The primary input is fed to a tapped delay line to generate tne adaptjve filter input signals. In thi s 

case, the input vector at the nth iteration becomes 

(for a filter of length N) 

and the weights are updated with the arrival of each new data sample 

x(n) 
Let xl (t) and X2(t), in Figure 2.3, be given as in equation 

(2.1) with the same assumptions. Let us consider the discrete time 

version of (1) for the following analysis. That is, 

and x2(n) = Si(n) + m2(n) = a Sl(n-D) + m2(n) 

where D is a positive integer less than N, representing the discrete 

time del ay between Sl (n) and S, (n) . 

The adaptive filter inputs are then given by 

x( n) = XI (n) = Sl (n) + ml (n) (2.26a) 

and d(n) = x2(n) =a Si(n-D) +m2(n) (2.26b) 

The adaptive filter, with inputs given by (2.26), can be thought of as a 

system attempting to insert a delay equal to the propagation delay between the two sensors in the primary filter input, x(n), a1 igning the 

signal component in time prior to subtraction to produce the error signal. Hence, one weight of the filter corresponding to the correct 

delay would be unity and all other weights zero. In practice, due to 

the fact that tne filter must interpolate between its discrete taps to 

provide delays that are noninteger multiples of the sample time, the 

weights converge to a shape that i s peaked at the correct delay. 
Sensor 1 

.Adaptive filter of length N I I 

Figure i.3. Modeling TOE using adaptive filter. 
Therefore, determination of the delay, as in the GCC method, requires 

estimation of the peak of the adaptive filter weight vector. 

To examine the shape of the weight vector, let us consider the 

frequency domain window corresponding to the least-mean-square set of 

weights, WLNS = ( vxx) -l Vdx . That i s, 

For d(n) and x(n) given by (2.26), this becomes 

Converting back to time domain and using the convolution theorem, the 

kth weight at the nth iteration is 

where 0 G K < N-1. Equation (2.28) is to be compared with equation 

(2.9), where W (K) is a discrete time version of R (T) obtained LMS X1 X2 

using the Roth processor. It is important at this point to note that 
equations (2.27) and (2.28) do not depend on G, (f) . Therefore, as 

2% 

suggested by Ahrned and Carter [17], it is desirable in this adaptive 

approach to have the signal with lower S/N ratio as the secondary input 

to the filter. From equation (2.27) we can see that when 

G (f) = const. x G (f) ml ml Sl Sl 

we have dLMS(f) = B e -jhfD9 where B is a real constant. In this 

case, equat i on (2.28) becomes 

Let us examine equation (2.30) for two types of input signals. 

First, consider the case when the input signals are broad-band, i .e. 

f f 

s s - - f < - where fs denotes the sampling rate. Evaluating (2.30) 

2 2 

for this type of signal yields 

= B sinc {T [K-D]) (2.31) 
where sinc (0) = sin(*) and 0 < K G N-1. In this case the weight 1. \ 

vector has the shape of a sinc function peaking with amplitude equal to 

one when K-D = 0 or K = D. The first zero crossing of (2.31) happens 

for 

Hence, the main lobe of the sinc function is only two resolutions wide, 

simul ating a delta function. However, as in the Roth processor, when 

the condition (2.29) on G,n (f) is not met, the main lobe of the sinc 

11 

function will be spread out by convolving with the integral part of 

(2.28). 

A second relevant case is one in which the input signals are band1 imited in frequency. Evaluating equation (2.30) for band-1 imited si gnals leads to a slightly different form of the weight function as discussed ~y Ahmed and Carter [17]. Consider an ideal band-limited signal 

as follows 

for wg<)wltwl 

GslS1(~) = 0 , elsewhere 

where the conversion to the radian frequency, w, was made for the ease 

of analysis only. For this type of signal, equation (2.30) becomes 
where 

= 48 ob sinc [ob(K-D)] cos [oc(K-D)] 

- @l *O - Wl*O Wb - - and oc - - 

2 2 

The peak in the weight function still occurs at the correct delay, that 

is, for K-D = 0 or when K=D, but the main lobe is spread out. That 

can be seen by examining the first zero crossing of equation (2.32). 

The first zero in the sinc function occurs for 

while the first zero in the cosine function occurs for 

Since illc 3 ub, the first zero in the cosine function occurs sooner 

'II- A than in the sinc and yields a main lobe width of 2 x - - -. So for 

2oc UJ C 

band-limited input signals, the spreading in the peak of the weight 

function is inversely proportional to the band width. For a white input 

A signal, a, = - and'tne main lobe width is two resolutions as pointed 

2 

out earlier. 
Other issues that should be considered in this approach include 

interpol ation, dynamic behavior of the filter, and the filter length. 

The time delay obtained from the estimated discrete weight function, 

WLMS(K), is an integer multiple of the sampling period while in practice the actual delay between the input signals may not be. To avoid 

the error introduced by the discrete delay estimate the weight function 

should be i nterpol ated to obtain its continuous counterpart WLMS(n ,K), , 

that is, 

Now, the desired noninteger sample of the estimated delay is given by 

the value of n at which WLMS(n ,K) is maximum. The study of the 

filter dynamic behavior by P. Feintuch, N. Bershad, and F. Reed [18], 

has shown a potential to track linearly changing time delays by just 

observing the peak weight move through the adaptive filter tapped delay 

line. It was also shown that the location of the peak of the weights 

1 ags the true time delay by an amount that depends on the delay rate, 

the signal correlation function, and the adaptive filter time response. 

Finally, the filter length should be chosen carefully to insure fast and 

yet accurate convergence. The lower bound on tne number of weights is 

depicted by the maximum propagation delay between sensors. From eyuation (2.32) we can see that as K increases the weight function 

W~~s (K) approaches zero. Thus, the filter length can be safely 
truncated in practice to the lower bound, in order to reduce the computati onal time. 

This section is concluded by a representative example i 11 ustrating 

the capability of the adaptive filter to detect the time delay after 

convergence of the weight vector. Let xi (n) be a stationary random 

process with a correlation function given by 

, 

R(n) = alnl with -1 < a < 1 

Assume that mi (n) = m2 (n) = 0 and x2 (n) = xi (n-0) where D is a 

positive integer given by 0 < D < N-1. The adaptive filter inputs are 

then, 

The input vector on the tapped-del ay 1 ine i s 

and the cross-covariance vector, Vdx, is 

I 

vdx = E [d(n)xi (n)] = E [d(n)~ (n) d(n)~ (n-1) *** d(n) x(n-N+1)] 
xT(n) = h(n) x (n-1) --* (n-N+l)] cb 

and the cross-covari ance vector, Vdx, would be 

The aut o-cov ar i ance matrix, Vxx, of ~(n) is 

and 
row, 

, 

i.e. the weight vector peaked at the 

correct del ay . 

2.6 Summary 

The general i zed cross correl ation approach was reviewed with more 

attention paid to tile Roth processor. The adaptive filter and its 

application to TDE problems was analyzed in detail. It was shown that 

the adaptive LIvISTDE method is equivalent to the Roth processor and tnat 

the adaptation process using the LMS algorithm, under stated conditions, 

ensures the convergence of the weight vector to the wiener-Hopf ~07ution 

that presents the best mapping between the two inputs of the filter. 

This mapping resulted in a weight vector peaking at the correct time 

delay with a peak resolution that is a function of the frequency band of 

the signals and their statistical characteristics. The advantages of 

the adaptive LMSTDE over the GCC method can be summarized as follows: 

- simp1 er imp1 ementation, since no FFT computat ions are i nvol ved 

- no prior knowledge of the statistics of the signal and noise are 

necessary 

- ability to track moving sources 
- 1 ess computational time. 

The advantages over the nonadaptive approach (i .e. direct implementation 

of the Wi ener-Hopf equation) are 

- simpler imp1 ementation, since no matrix operations are involved 

- the ability to track moving sources 

- the filter length can be easily increased in order to decrease 

tne bias in the estimate. 
Chapter 3 

SYSTEM DESIGN 

3.1 Introduction 

The theory behind t ime-del ay estimation techniques has been revi ewed in Chapter 2. This theory can now be used in the design process of 

the system. Based on the anticipated input signals and assumptions that 

can be made about the statistical characteristics, the choice of a suitable TDE technique is made in section 2. Section 3 then presents tne 

implementation of this technique and relevant design issues in tne 

specific application on nand. The system specifications and the way 

they are related, derived, or determined are the subject of section 4. 

Finally, the design process is summarized in section 5. 

3.2 Design Assumptions 

The choice of the TDE approach and the assumptions on which it is 

based are the subject of this section. Not much study has been dedicated to the high frequency band (I to 20 Hz) of infrasounds in the 

past, especi a1 ly for sounds re1 ated to meteorological events. Thus, 

little is known about the statistics of these infrasounds and assumptions have to be made with regard to recent investigator's predictions. 

These assumptions will be listed below with a relevant discussion about 

each of tnem. 

a. Single Source. More than one meteorological event can take 

pl ace at one time producing a mu1 tiple delay between sensor 
outputs. The potential to separate delays and identify the 

different sources is feasible but it is outside the scope of 

this study. So the design process is restricted to a single 

source of sound waves. 

b. Stationary Si gnal . The stationary assumption i s safely made 

because the data file spans a time window that is only a fraction of the even$ Is 1 ifetime. The l ifetime of a meteorological 

event can be anywhere from a few minutes, in the case of a 

microburst for example, to several hours in the case of a frontal system [I], [19]. Tnis time frame is to be compared with 

the file time that, as will be seen later, is 20 seconds. 

c. Lhcorrel ated Signal and Noise. Infrasounds emitted by a re1 atively remote source cannot be expected to correlate with local 

noise present at the microphone array. 

d. Uncorrel ated Noise at the Different Sensors. This assumption 

can be justified by having a long array foot so that local 

noise at one sensor is not swept across the array. 

e. Frequency Band. The low frequency band (10-6 to about 1 Hz) of 

the infrasound spectrun is more related to static pressure 

fluctuations that are part of a turbulent motion. Since in 

this project, the intenti on is to detect acoustic (propogatiny 

wave) energy and as 1 ittle as possible of the turbulent pressure fie1 d, the frequency band should be in the nigher part of 

the spectrum. The 1-to-16-Hz band has been experiment a1 1 y 

studied by Eric S. Posmetier [20]. He concluded that "the 

1-to-16-Hz i nfrasound may be radiated by clear air turbulence 
and may be a basis for a remote passive detection system." 

Based on the physics of i nfrasounds and Posmetier's prediction, 

a signal frequency band of 1-to-20-Hz is anticipated. 

f. Spectral Shape of Signal and Noise. Based on past experience, 

the signal and noise spectra are expected to peak in the range 

.5-to-3 Hz. 

Now that the assumptions about the data have been established, the 

choice of a suitable time-delay estimation approach can be made. The 

data statistical characteristics required for the General i zed CrossCorrelation method have been assumed in a, b, c, and d. Thus, the GCC 

method is applicable in this passive detection problem. As was shown in 0 

section 2.2, the basic cross-correl ation has an undesirable spreading in 

its peak, which can be improved by a frequency weighting. The weighting 

proposed by Roth leads to a weighted correlation given ~y 

For the assumption made in f, about a similar spectra1 shape for signal and noise [i .e., Gm,(f) E const. GSS(f)], this correlation 

becomes 

Thus, the Roth processor is appropriate for the assumed data, since its 

output is a delta function with a magnitude equal to one at the correct 
delay. Implementing the Roth processor in the frequency domain, however, has the disadvantages discussed in section 2.3. & the other 

hand, the advantages discussed in section 2.6 make the choice of implementing this processor in the time domain, using the LMS adaptive algorithm, evident. The choice of the TDE approach has been made. The 

question now is whether the adaptive filter discussed in section 2.5 is 

totally suitable for this application. This is the subject of the next , 

sect i on. 

3.3 Filter Desiyn 

The filter structure discussed by F. Reed, P. Feintuch, and N. 

Bershad [16], and presented in section 2.5, requires prior lead-1 ag 

information about the monitored signals in order to successful ly est imate the delay. The delay estimate can be made only if the leading 

signal is fed as the primary input to the filter. This can be seen by 

examining Figure 2.3 where the input vector 

component of the second sensor output, d(n) must 1 ag x (n) so that 

the delay estimate can be made. 

In acoust i c detection problems, the 1 ead- 1 ag i nformat i on i s 

unknown. T'nis fact results in the necessity of applying the algorithtn 

twice to the same data, swapping channels the second time, in order to 

get the lead-lay information as we1 1 as the delay estimate. This 

11 difficulty can be solved by introducing an - point delay in the 

2 
secondary f i 1 ter input as i 11 ustrated i n Fi gure 3.1. Introducing thi s 

delay in d(n) leaves enough leading as well as lagging information in 

X(n) to be compared with d(n) to eventually produce the estimate. 

, In this implementation, the filter inputs will have a slightly 

different form than before. With discrete time signals given as in 

Chapter 2, 

, 

X2 (n) = S2 (n) + rn2 (n) = a Sl (n-D) + m2(n) 

the filter inputs are now given by 

and d(n) = ~2 (n-N') = a Si (n-N' -D) + m2(n-N') (3.lb) 

N where 0 G n G record 1 ength, I Dl c N' , N' = integer (-) , and M i s an 

2 

odd positive integer equal to the filter length. If, as assumed 

earlier, 9 (n) and m2(n) are uncorrelated, having %(n) delayed by 

N' points have no bearing on the analysis. A1 1 the analysis, presented 

in Chapter 2, of the weight vector is still applicable here except that 

D must be replaced by D+N1 everywhere in the derivations. Equation 

(2.31) becomes in this case 
Sensor 1 1 

1 

Sensor 2 

I Adaptive filter of length N 

* N=integer(y) 1- - -- - - -- - -- - - - - - - - - - - - - 

2 

Figure 3.1 TDE lead-1 ag model . 
w' (K) = 6 sinc{r [K - (u+N)]} LMS 

As can be seen from (3.2), the weight vector now peaks at K = I)+N'. So 

when 0=0, the weights peak at N' (the middle weight), and the peak 

is shifted right or left with positive or negative values of the delay, 

respectively. 

In section 2.4, it was shown that convergence of the weight vector , 

to the Weiner-Hopf solution can be obtained only if the stability constant Ks is bounded by equation (2.25), that is, 

Input power 

However, mai ntai ni ng a constant Ks through the adapt at ion process can 

lead to instability in situations where a rapid cnange in the power 

level of the signal may occur. For this reason, Ks is updated in this 

design with every iteration of the algorithm. An estimate of tne input 

power of the signal, based on the signal cornponents in the tapped de1 ay 

line, is made with the arrival of each data sample. The value lKsJ is 

then chosen to be less than the inverse of the power estimate by a factor to be experimentally determined. Updating Ks in this fashion ensured a fast, stable convergence of the weight vector regardless of the 

power level of the input. 

The final design issue is the convergence criteria. Due to tne 

fact that our sampled data records may not contain any coherent signal s 
all the time, and due to the fact that the shape of the weight vector is 

unpredictable in these situations, a convergence criterion must be derived to distinguish between coherent and noncoherent input data before 

attempting to estimate the delay. This criterion is proposed as 

follows. Let us monitor the peak of the weight vector at every iteration as the data record is processed through the algorithm. By storing 

the abscissa value (corresponding to the discrete time del ay) of the 

peak with every iteration, a data sequence i s produced that can be plotted to produce a TDE vs. time plot. A variance criterion can then be 

applied to this TDE vs. time plot to determine whether a convergence 

has occurred or not. The data sequence representing tne abscissa values 

of the peaks is as long as the data record itself. Thus, for snort-term 

stationary data, the variance criterion should be applied to a sliding 

window, of lengtn less than the data records, through the sequence instead of computing the variance for the whole sequence at once. Whenever the variance criterion is met, the delay estimate is made based on 

the abscissa values contained in the window where the convergence has 

occurred. This provides a better estimate of the delay, compared to 

taking the abscissa value of the peak in the last iteration, especially 

in noisy environments. 

A detailed description of our implementation of the adaptive LMSTDE 

algorithm is presented next in Figure 3.2, Figure 3.3, and Figure 3.4, 

by means of flowcharts. In Figure 3.2, it is assumed tnat the sampled 

data records nave been stored in arrays x and d of length L, wnere 

L is the data record length. W is an N point array, where N is 

the filter lengtn, storing the weight vector components. The abscissa 
0 Start 

Update p as II = v + x2(n) - x2(n-~+1) 

.1 

-1 Compute Ks = 

N-1 

Initialize: n = 0, ,, = z x2(i), w(i) = 0 for 0 5 i 5 N-1 . i =O 

" I N-1 Compute sum = z X(n-i) w(i) 

i -0 

<, 

Compute E = sum - d(n-N)) 

- - \1 

Update the weights as w(i) = w(i) - 2 K~ x(n-i), 0 i ~-1 

+ 

Apply peak detection routine & store abscissa of peak in p(n) 

1 

Compute variance 

I 

1 Compute delay 

I I 

Figure 3.2. Li4STDE algorithm's flowchart. 
0 start 

Least = 1000 , pt = 1000 

I 4 

I Cmpute ,,I = + P(i), p2 = u2 + p2(i) for 0 s i K-1 

I 

Least = VR, Pt = i I 

VR = 'Least I_I 

Figure 3.3. Variance computation. 
Start Q 

Compute D = D + P(i) for (Pt-K i s Pt-1) 

h I I 

Figure 3.4. Delay estimate computation. 
values of tne peaks will be stored in an L point array P. VK and 

VT denotes the vari ance and the vari ance thresh01 d, respectively. The 

estimate of the stability constant, Ks, is computed with every iteration as 

where C is a predetermined factor and 

N-1 N- 1 

u = 'Z ~[x:(n)] = N x E [X2(n)] = I X2 (n-i) 

i =O i =O 

Tine peak detection block in Figure 3.2 is a simple search routine, wnil e 

the variance computation is described in Figure 3.3 and tne delay cornputation in Figure 3.4. The variance computation is made, as described 

earlier, over a sliding window of length K, where K < L. The equation used for the variance estimate is 

Finally, the delay estimate is made by simply averaging the delay values 

contained in the window that yielded the minimum variance level. 

By inspecting the algorithm flowcharts, it is easily seen that the 

computational complexity is proportional to L x N. In practice L is 

the minimum number of iterations necessary for the convergence of the 

weight vector to the diener-tiopf solution. Therefore, for a specific 
statistical type of data, L is a constant and the complexity becomes 

of order N. This design, accommodating positive and negative delays, 

impl ies a lower-bound on N of twice the maximum propagation delay 

between two sensors. Compared to the filter discussed in section 2.5, 

the filter length now is twice what it was before. With the complexity 

being of order N, it appears as if no speed improvement has been 

established over the previous design, where two passes were needed to 

obtain the de1 ay estimate. Actual impl ementati ons, however, proved 

about 10% time saving. The time saving was mainly obtained oy avoiding 

the overhead computation i nvol vea in every pass as can De seen in the 

flowcharts. 

3.4 Specifications and Hardware 

This section presents and justifies the specific parameter values 

in the algorithm implementation as well as in the hardware system used. 

These specifications are summarized in Table 3.1. They can be divided 

into three categories. The first category includes parameter values 

derived from the assumption made in the previous section about the signal frequency band. Tne second consists of parameters that were experimentally determined, while the ones related to hardware make up the 

third category. 

The starting point in tne first set of parameters is tne frequency 

band tnat was specified as 1 to 20 Hz. An immediate result of this band 

is the sampling rate. A sanpling rate of 40 iiz, satisfying the ~Vyquist 

rate, will not be adequate, tne reason being that the signal has to be 
Table 3.1. Parameter Specifications 

low-pass filtered at 20 Hz with a non-ideal rol 1-off filter. To el imi - 

nate any alias effects produced by attenuated frequencies move 20 dz, a 

sampling rate of 51.2 Hz was cnosen. The exact value of 51.2 Hz is 

desirable because it produces a rounded frequency resolution of .2 Hz on 

a 256-point FFT spectrum. Given the sampling rate, the array foot 

1 ength (distance between two sensors) becomes a compromi se between 

computational speed and beari ng resol uti on. The shorter the array foot, 

the less the (maximum) delay, and the lower the resolution. The filter 

length, and hence the speed is directly proportional to the maximum 

propagation delay as was seen in the previous section. The bearing 

VALUE 

2 Hz - 20 Hz 

51.2 rlz 

800 ft. 

12 bits 

2.4 mV/bit 

1024 points 

128 points 

77 points 

10 

1 

48 

PARAlYETER 

Frequency Band 

Sample Rate 

Array Foot 

A/D 

, Quant. Level 

Fi le Length 

Record Length 

Fi 1 ter Length 

Scale of Feedback 

Estimate 

Variance Threshold 

Vari ance Window 

IJuTATI ON 

- 

s 

R 

- 

- 

M 

L 

N 

C 

V T 

K 
resolution is related to the array foot as 

A9 = ~o"-cos-~ ( At x speed of sound 

array foot length 

1 

where A0 is the bearing resolution and At is the sampling period. 

To develop insight for the range of values involved, we can compute the 

filter length and bearing resolution for three different values of the 

array foot. (ht = 1/51.2 sec, sound speed = 1080 ft/sec, filter 

length = 2 Omar+l) 

Array foot 2000 ft. 800 ft. 100 ft . Filter length 200 77 10 

deari ng resol ut i on .6" 1.5" 12" 

Tie first set of values implies a good resolution at the expense of slow 

speed while the third set implies a high speed computation at the expense of poor bearing resolution. The compromise is obtained in the 

middle set. Note that the chosen distance between sensors is not likely 

to violate our assumption of uncorrelated noise across the array. This 

is true since infrasound noise produced mostly by local wind fluctuations and small eddies are not expected to propagate over an 800 foot 

distance. 

Experimental 1 y determined parameters are the record 1 ength (nunber 

of iterations), scale of stability constant, variance window, and variance threshold. Tie first two parameters, and the second two, are related. The record length, L, should be the minimum number of iterations that leads to convergence in the weight vector to the Niener-Hopf 
solution. The rate of convergence is control led by the stability constant, KS, as was seen in equation (2.21) in section 2.4. Keeping one 

of the two constant and varying the other, the system performance was 

observed and the values were chosen as 

L = 128 iterations 

, C = scale of Ks = 10 

The system performance, as L and C vary, is simulated in section 

4.2. The variance window, K, and the variance threshold, VT, were 

used in the convergence criteria to determine whether the data is coherent or not. Ideal response would have a zero variance across a 128- 

point (record length) window, but obviously that is too tight a criterion in practice. Based on the response of the algorithm to some experimental data gathered, incl uding space-shuttle 1 aunch noise, moderate 

values of 48 points and 1 were chosen for K and VT, respectively. 

It is worth mentioning at this point, that although the algorithm 

computation is Dased on 128 iterations, the data files sampled are 1024 

points each. Tnis file length is necessary for subsequent spectral 

anal ysi s . 

A block diagram of the system hardware and most of its specifications is presented in Figure 3.5. The type of sensors used in the array 

is the GLOBE 100 C Capacitor Microphone by GUS Manufacturing Inc. The 

high sensitivity ( .2 vol t/dyne/cm2) and the low frequency response ( .1 to 

(.I to 500 Hz) on these Microphones make them very adequate for this 

application. Their dynamic range is 74 dB, .004 to 20 volts peak-to-peak 
1 GLOBE 

Microphone 

Array 

IB&K 2 Hz, 2 pole 20 Hz, 16 pole 12-Bit ' 

Bessel Computer Butterworth A/ D 

Low-Pass and 

Amp1 ifier High-Pass 

Filters Converter Peripherals 

, , , , 

, Filters 

, , , i/ 

Attenuator 

Figure 3.5. System block diagram. 
peak. This dynamic range implies an electrical noise amp1 itude of 

(upto) 4 m volts p-to-p. To avoid any electrical noise in the sampled 

data, the analog-to-digital converter should have a quanti zation 1 eve1 

of at least 2 mV/bit when no amplification is used. But since a large 

dynamic range A/D is obviously desirable, a 12-bit, 2.4 mV/bit AID converter was chosen. An analog input range of -5V to +5V corresponds to a 

digital output range of 0 to 4095. Therefore, as long as no ampl ification is used on a microphone's output, the data are free from electrical noise. Tie specifications on the remaining nardware blocks are 

totally independent. The array configuration was intended to provide 

sy~metry between the microphones so that delay-to-angle conversion can 

be made easier. Due to the experimental nature of tnis project and due 

to the lack of knowledge about the sound pressure level of infrasounds 

generated by meteorological events, the ampl ifier/attenuator block was 

necessary. Even though a signal frequency band of 1 to 20 Hz is anticipated, the input is high-pass filtered at 2 Hz in order to eliminate the 

strong wind noise in the region of .1 to 2 Hz. Finally, the choice of 

the APPLE computer, as the system processing unit, was made to meet the 

economi cal , port abl e system requirement. 

3.5 Summary 

The design process started by evaluating the anticipated signal and 

noise. The main assumption made was that of similar spectral shape for 

signal and noise over a frequency band of 1-to-20-Hz. This assumption 

led to the choice of the Roth processor which, in this case, produced a 

snarp peak at the correct delay in the weighted correlation function. 
Examining the facts of section 2.6, proved that implementing the Roth 

processor in the time domain (rather than frequency domain as in the GCC 

method) using the LMS adaptive filter is advantageous. The specific 

N filter design for this application was then presented. An - -point 

2 

del ay in the secondary filter input was necessary to accotnmodate for 

positive and negative delays. Updating tne stability constant, I(~, 

with every iteratiaon was made for fast sta~le convergence of tne weignt 

vector when the power level of the signal may change rapidly. To discriminate against non-coherent input signal s in the automatic detection 

,node of the algorithm, a convergence criterion was devised. Tnis criterion was based on the variance level of the delay estimate with time. 

The design phase ended by schematically presenting the exact implementation of the design and specifying the software, as well as the hardware, 

parameters values. 
Chapter 4 

SIMULATION AND RESULTS . 

4.1 Introduction 

Aided by the theory of TDE techniques reviewed in Chapter 2, the 

system design has been described in the previous chapter. This chapter 

presents an evaluation of the foregoing design. The response of the TOE 

algorithm, as presented in section 3.3 is evaluated in section 2. Its 

response for determi ni sti c data sequences is simul ated as the different 

parameters are individually varied. Its effectiveness in tracking time 

delays between real, long-range, infrasounds is also illustrated. kction 3 describes the complete software package, with the TDE algorithm 

being its central processing element, used in the actual field operation 

of the system. Its routine operation and genera1 results over the 

period of June to September 1984, are also discussed. The system is 

evaluated with respect to man-made and weather-related infrasounds. The 

man-made infrasounds are discussed in section 4, with an illustrative 

example using signals from jet-powered aircrafts, while the weatherrelated infrasounds are the subject of section 5. 

4.2 System Behavior 

This section is dedicated to the study of the system response for 

deterministic input data as well as real infrasound waves generated oy a 

far-field source. Deterministic data are used to study the system 
Denavior as the different parameters, with effect on the response, are 

varied one at a time. Tne type of data used in the simulation and tne 

response for a nominal set of parameters are discussed in sub-section 1. 

The set of parameters can be divided into two groups. The first group 

consists of the adaptive filter parameters that were experimentally . 

determined in the design. This group is the subject of subsection 2. 

The second group, contaiqing the signal parameters (signal characteristics), is discussed in subsection 3. Finally the response for real 

signals consisting of the space shuttle launch noise is illustrated in 

the fourth subsection. 

4.2.1 Nomi na1 Response 

The data used in this simulation are white data sequences generated 

on an Apple I1 computer using a random nunber generator function. Nhite 

data are used for both signal and noise. Even though data with decaying 

spectra seem more appropriate for tne simulation, white sequences were 

chosen for the ease in which they are generated. As explained in section 2.5, both types of data are expected to produce tne same response 

since the weight vector, for signal and noise having the same spectral 

shape, is a sinc function regardless of what the spectral shape is. A 

256-poi nt sequence and its corresponding auto-power spectrum, representing the type of white sequences used, are plotted in Figure 4.1. The 

auto-power estimate is obtained by a 256-point FFT, and represented by 

the first 129 components since the 1 ast 127 are a mirror image of the 

first. The spectrum is expected to be completely white, in contrast to 

the one in the figure, as the nmber of overaged 256-point periodograms 
I I I 1 I 1 I I 

0 64 128 192 255 

Sample Number 

(a 1 

Di screte Frequency 

Figure 4.1. Type of random data sequences used for simulation. 

(a) 256-point sequence. (b) Its auto-power spectrum. 

\ 
approaches i nf i ni ty . 

The nominal response of the algorithm for positive and negative 

delays is presented in Figure 4.2 and Figure 4.3, respectively. The 

positive delay in Figure 4.2 was generated in the input sequences as 

fol 1 ows: 

where x (n) is the primary input to the filter, d(n) is the secondary 

input, Rl (n) , R2 (n) and R3 (n) are uncorrel ated random sequences. The 

embedded delay is four samples while the other parameters are set to 

their nominal values as specified in the figures. The negative delay 

was obtained by simply swapping the adaptive filter inputs, i .e. : 

Tne fact that ~(n) now 1 ags d(n) by four samples leads to a negative 

delay in the output of the filter as illustrated in Figure 4.3. 

Tine first plot in both figures was obtained by monitoring the 

abscissa value of the peak in the weight vector at every iteration of 

the algorithm. The second plot is the final set of weights. The 
Nuniber of iterations 

Weight number 

Figure 4.2. Nominal filter response for D = +4 samples. (C=10, 

N = 17, L = 128, SilR = 1). (a) TDE vs. number of 

iterations. (b) Final weight vector. 
Number of iterations 

Weight number 

(b) 

Figure 4.3. Nominal filter response for D = -4 samples. (C=10, 

N.17, L.128, SNR=l). (a) TDE vs. number of iterations. 

(b) Final weight vector. 
capability of tracking positive and negative delays is evident from the 

figures. The "steady-state" delay has shifted from +4 to -4 in the TOE 

vs. time plot, while the peak in the weight vector has shifted right and 

tnen left by 4 units as expected. These plots are to be compared with 

the response obtained in the next two subsections where the various 

p aramet er s are v ar i ed . 

, 4.2.2 Fi 1 ter Parameters 

The filter parameters include the filter length (N), stability 

constant scale factor (C), and the record length (L) . In the design 

process, a lower bound on N was set while C and L were experimentally determined. This simulation verifies the necessity of tne lower 

bound on N and explores the basis on which C and L were determined. The input sequences used in this subsection are the ones previously given in section (4.1). 

The lower bound on tne filter length (number of weights) was derived in section 3.3 as (2xDmax+1). Twice the niaximun delay was necessary since the filter accommodates positive and negative delays. A 

number of weights less than the lower bound will obviously give no estimate. A very large nunber of weights, on the other hand, will not 

improve the response since the sinc function (the shape of the weight 

vector) approaches zero as the nmber of weights is increased. This 

behavior is exemplified for D=4 in Figure 4.4 with N=7 and in Figure 

4.5 with N=33. The remaining parameters are set to their nominal 

values in both figures. The response in both cases is to be compared 
0 3 2 64 96 128 

' Number of iterations 

Figure 4.4. H=7. Delay can not be estimated with Nc2xDmax+l. 

(D=4, C=10, L=128, SNR=l). (a) TDE vs. 

number of iterations. (b) Final weight vector. 

+, 

- 

- 

- 

- 

- 

- 

- 

- 

0 I 

- 

- 

- 

- 

- 

- 

- 

L~~III~I~I~III~~~ 1 IIIIIIIitIirfiIiL 

-64 -32 -16 0 16 32 64 

Weight number 
0 32 64 96 128 

Number of iterations 

(a 

Weight number 

Figure 4.5. N=33. Ho improvement over N=17. (D=4, C.10, L=128, 

SNR=l). (a) TDE vs. number of iterations. (b) Final 

weight vector. 
with the nominal response in Figure 4.2 where N=17. 

The effect of varying C on the stability of the system is demonstrated in the next two figures. It was seen in section 3.3, the bo~nds 

on KS are given by: 

The value used in this implementation is given as 

where C is a positive scale factor that has to be greater than. one to 

satisfy the stability condition. The output of the filter is plotted in 

Figure 4.6 and Figure 4.7 for C=.5 and C=200, respectively. Having 

C less than one led to an unstable system as predicted by Widrow [ll], 

while increasing its value behind ten hardly affected the response. 

The record length (nunber of iterations) was chosen as 128. The 

purpose of this simulation is to show that a 128-iteration process is 

sufficient for the convergence of the weight vector to the Wiener-Hopf 

solution. This is illustrated in the next two figures in comparison 

with Figure 4.2. Figure 4.8 illustrates the ambiguity in the filter 

output caused by too small a number of iterations. The expected 
' Number of iterations 

Weight number 

Figure 4.6. Unstability for C-0.5. (0=4, N-17, L-128, SNR-1). (a) TDE vs. number of iterations. (b) Final weight 

vector. 
' Number of iterations 

(a> 

Weight number 

Figure 4.7. C.200. No improvement over C=lO. (D=4, iY.17, L=128, 

SNR=l). (a) TDE vs. nunber of iterations. (b) Final 

weight vector. 
' Number of iterations 

(a> 

Weight number 

(b 

Figure 4.8. L-32. Ambiguity in filter output due to short record 

length. (D=4, N=17, SNR=l) . (a) TDE vs. number of 

iterations. (b) Final weight vector. 
response, after a satisfactory convergence, was obtained for L=128 in 

Figure 4.2. A larger number of iterations produced sligntly improved 

results, though not appreciable, for the additional computation. This 

aspect is shown in Figure 4.9 for L=512 points. 

4.2.3 Signal Parameters 

This subsection illustrates the system response as the cnaracteristics of the filter input; are changed. The characteristics considered 

here are signal-to-noi se power ratio (SNH), signal frequency bandwidtn, 

and time varying delays. 

The signal to noise ratio used in the nominal response of Figure 

4.2 was SNR=l. The response is evaluated in the next two figures for 

SNR =.5 and SNR=2. The input sequences were generated as 

x(n) = SxRi (n) + R2 (n) 

where x(n) is the primary input to the filter and d(n) is the secondary one. Rl (n) , R2 (n) , and R3 (n) are uncorrel ated random sequences 

with the same amplitude. The scale factor S equals .5 in Figure 4.10 

and 2 in Figure 4.11. In all three cases, the response of the system 

was satisfactory after 128 iterations. However, as can be seen from 

the figures, the lower the SNR the slower the response (convergence). 

The next simulation demonstrates the dependency of the main lobe 
0 32 64 96 128 . Number of iterations 

(a> 

Weight number 

(b 

Figure 4.9. L=512. Improved response over L=128 but not necessary. 

(1)=4, C=10, N=17, SNR=l). (a) TDE vs. nunber of iterations. 

(b) Final weight vector. 
Number of iterations 

(a> 

W'eight number 

Figure 4.10. SNK=1/2. Slow Convergence. (0=4, C=10, t4=17, L=128). 

(a) TDE vs. nunber of iterations. (b) Final weight vector. 
Number of iterations 

(a 1 

Weight number 

Figure 4.11. SNR.2. Faster convergence. (0.4, C-10, N-17, L-128). 

(a) TDE vs. number of iterations. (b) Final weight vector. 
width of the weight vector on the signal frequency band. The expression 

for the weight vector for an ideal band-limited signal and white additive noise was derived by Ahmed and Carter [17] as 

WLMS(K) = 4Bwb sinc[ub(K-D)] cos [wc(K-d)] 

where = Wl-00 6 is a real constant, wb - , w = Wl-wo 

C W1 and wo 

2 2 

are the signal frequency 1 imits, i .e. Gs s (w) = 0 for w <wo and w 

> w1 , where wl > Wo . It was shown in section 2.5 that the width of 

IT ZIT the main lobe of the weight vector is given by -, which becomes - 

for wo = 0. The input sequences were generated as 

x(n) = S x Hl (n) + Re(n) 

where (n), d(n), K2 (n), and R3 (n) are the same as before, S is a 

scale factor used to obtain a SNR=1, and R (n) is an ideally filtered 

sequence specified by its amp1 itude spectrun as 

In Figure 4.2, the frequency band of Kl (n) was given by wo = 0 and 

w, = IT , i .e. white sequence. The main lobe in that case was 2 
& resolution ( = 2) wide. For the simulation in Figure 4.12, the 

w1 

71 bandwidth of Rl(n) is limited to -, i.e. 

4 

n 71 

G~l ~1 

(w) = 0 for - <w c n + wo =O and wl =- 

4 4 ' 

2n As seen in the figure the main lobe width now is - - 8 resolutions. 

, - IT 

4 

n Finally, wg = 0 and wl = - were used in Figure 4.13. Producing a 

8 

main lobe that is 16 resolutions wide. 

It was mentioned in section 2.6 tnat one of the advantages of the 

LMSTDE algorithm over other methods is that it has the aoility to trac~ 

moving sources. The routine is tested here for this ability by linearly 

varying the delay between the input signals in a noisy environment. The 

input sequences were generated as 

where Rl (n), R2 (n) and R3(n) are uncorrelated random sequences as 

before, and D(n) is a time varying delay. The signal-to-noise ratio 

used here is 2. In Figure 4.14, with U(0) = 0, D(n) was increased by 

one every 64 points in 896-point input sequences. The solid line 
I 

Number of iterations 

(a) 

Weight number 

(b) , 

IT Figure 4.12. W1= - . Main lobe width of 8 resolutions. (0.4, C=10, 4 N.17, L.128, SNR=l), (a) TOE vs. nunber of iterations. 

(b) Final weight vector. 
Number of iterations 

(a) 

Weight number 

(b I 

n Figure 4.13. Wl= . Main lobe width of 16 resolutions. (0~14, C=10, 

N=33, L=128, SNR.1). (a) TDE vs. nunber of iterations. 

(b) Final weight vector. 
(suo~qn LosaJ u~) Ae lap aqa~3s La 
represents the actual delay embedded, while the dashed line represents 

the instantaneous estimate of the delay as the abscissa value of the 

peak in the weight vector during 896 iterations. Figure 4.15 differs 

from Figure 4.14 in that D(n) is varied at the rate of 1 point every 

128 points instead of 64. The adaptive filter tracked more closely the 

actual delay in the second case. This is to be expected, since the 

weights were given more .time to converge. In both cases however, D(n) 

represents a relatively fast varying delay simulating a fast moving 

source. Considering this fact, one can conclude that the system 

response is acceptable and that the time 1 ag, predicted by Feintuch, 

Bershad, and Reed [18], between the estimate and the actual delay is not 

as severe as it appears. 

4.2.4 Shuttle Launch Noise 

Long-range infrasound signal s produced by a 1 arge rocket 1 aunc h 

have been recorded and studied by several investigators [22]. These 

signals are similar in their statistical characteristics to weatherre1 ated infrasounds. Thus, they can be used to simul ate the anticipated 

infrasounds generated by a f ar-f ield meteor01 ogical event. In thi s 

subsection, the adaptive filter algorithm is applied to the space 

shuttle launch signals and the overall response is evaluated. 

Space Shuttle VIII was launched from Cape Canaveral, Florida, at 

2:31 a.m. EST on 30 August 1983. About 52 minutes later the launch 

signal arrived at Wallops Island, Virginia, and was sampled on an array 

of four infrasonic microphones at the rate of 14 Hz. The array configuration is illustrated in Figure 4.16. Time histories and auto-power 

spectrums, of one channel before, during, and after the event, are presented in Figure 4.17 and Figure 4.18, respectively. These two figures 

represent a verification of the event and show the background of the 

signal that can be considered as the uncorrelated noise on the different 

microphones. The adapti ye f'iJ ter algorithm was appl ied to three pairs 

of channels during a 9.14 sec. (128 points) time frame. The output of 

the microphones during this time frame is plotted in Figure 4.19. The 

adaptive filter output for C2 x C3 (channel 2 cross channel 3), C3 x C4, 

and C4 x C2 is shown in Figure 4.20, Figure 4.21, and Figure 4.22, 

respectively. Making the estimate of the delays based on a 48-point, 

zero variance, window in the TDE vs. time plots yields the following 

results: 

C2 leads C3 by .455 sec. (+7 resolutions) 

C3 1 ags C4 by .780 sec. (-12 resolutions) 

C4 leads C2 by ,325 sec. (+5 resolutions) 

dased on these delays and on the array configuration, a source direction 

of about 210" from north can easily be obtained using trigonometric 

identities. This estimate of the bearing is to be compared with the 

true direction of 207" as illustrated in Figure 4.23. The error in the 
1 

I ; 

/ 

750 ft. I 

4) 2 

1 

I 

/ 1 / I 

Figure 4.16. Array configuration used in sampling STS-8 launch 

signal . 
1 I I I 1 1 I I 

9.14 seconds 

Discrete time 

I 1 I t I I I I 

9.14 seconds 

.r 

C 

Oi sc,rete time 

- 

I I I I I I t 

Figure 4.17. Time histories of STS-8 launch signal (in equivalent 

arbitrary units) on one channel, (a) before, (b) during, 

and (c) after the arrival of 1 aunch signal. 

9.14 seconds 

Discrete time 
0 

lllllll~ll~~~~llllltlll~lt~ IIIII 

0 7 Hz Di scre-te frequency 

Discrete frequency 

Discrete frequency 

Figure 4.18. Auto-power spectrums of STS-8 1 aunch signal (in arbitrary 

units) on one channel, (a) before, (b) during, and (c) 

after the arrival of launch signal. 
I I t I 1 1 I I 

9.14 seconds 3i screte time 

I I I I 1 I 1 I 

9.14 seconds 

Discrete time 

a, 

-0 

3 

C, 

.r 

Figure 4.19. Time histories of STS-8 launch signal (in equivalent 

arbitrary units) on three inicrophones, (a) channel 2, (b) 

channel 3, (c) channel 4. 

- . . . . -. ;: . . ............ . . . . ... ; '., . :?. -* .-*. ., . . : , ; ,:.*p,.iv.. ..... .-.,: . : .. .., % ,,..*;'.-,. .... -3 ' ... ' ..... g .: - ,. . . . . .- .. .. 

I I I I I I t 

9.14 seconds 

Di screte time 
, Number of iterations 

Weight number 

(b 

Figure 4.20. Filter response for launch signal, C2 x C3, (a) TDE 

vs. time, (b) Final weight vector. 
itumber of iterations 

(a 

Weight number 

(b) 

Figure 4.21. Filter response for launch signal, C3 x C4, (a) TDE 

vs. time, (b) Final weight vector. 
 umber of iterations 

(a) 

Weight nur~iber 

Figure 4.22. Filter response for launch signal, C4 x C2, (a) TOE 

vs. time, (b) Final weight vector. 
Figure 4.23. True direction of shuttle launch site. 
estimate is within the bearing resolution. In this case, the bearing r I 

resolution (for fs = 14 Hz and array foot of 650 ft) is about 6.8'. 

4.3 System Operation and Results 

The previous section demonstrated the effectiveness of the LMS 

adaptive filter algorithm in estimating the delay between the outputs of 

two sensors. The algorithm, as discussed in section 3.3, can now be 

used as a subroutine in d 1 arger software package for a complete, actual 

field operation of the system described in section 3.4. This section 

presents a description of this software package, the mode in which it is 

operated, and the way the results are evaluated. 

The software package is depicted in Figure 4.24. It consists of an 

automatic detection routine, as we1 1 as a bearing computation, of infrasonic sound waves. The whole process is automatic in the sense that no 

user intervention is needed after the program is started and the pdrameters are set. However, this automation is 1 imited by the data storage 

capacity of the computer. The process is started by simultaneously 

sampling three 1024-point files from the three sensors* in the array. 

The data files are stored in temporary buffers until they are analyzed 

by the LMSTOE algorithm. The algorithm, as described in section 3.1, is 

then applied to each pair of data files at a time, considering only the 

first 128 points as the filter data record. If the convergence test, 

using the variance criterion, fails for any pair of sensor outputs, the 

process will be stopped. and a new set of data files is sampled. When 

*Only three out of the four sensors in the array are at this stage. 
0 Start 

L 

I Set parameters 

m Sample three channels 

\Ir 

Apply adaptive filter to all three pairs of channels 

I 

Store data LrJ 

' Evaluate angles 

P Print delays and angles 

Figure 4.24. Software flowchart. 
all three pairs of files show good coherent data through acceptable 

variance levels, the data files are transferred to a floppy disk for 

permanent storage. 

As the detection phase ends, the spectral content of the signal, 

the bearing, and the nature of the source are then evaluated. The 

analysis begins by testing the delays, obtained from the LIASTDE algori thm, for consistency assuming plane wave propagation and 0" elevation. 

If the consistency test is not passed the delays are printed for later 

analysis; otherwise three angles, corresponding to the three delays, are 

evaluated using simple trigonometric relations. The delays and angles 

are then printed for the record, and the whole process is started over 

again. 

The consistency test discriminates between plane waves and cylindrical, or spherical propagating waves. The 1 ater category may occur if 

the source is in the near field of the array (within one mile from the 

center of the array). Since the intention of this research is the early 

detection of far-field events, this category is neglected by the software. The delays obtained in such cases are only printed for direct 

analysis by the user. The trigonometric relations used to evaluate the 

bearing in case of a plane wave were formulated for hori zontal ly incoming waves only. In this case, the three computed angles will be the 

sane while they will differ from each other for a plane wave propagating 

with a non-zero elevation. The difference among the angles can be 

easi-ly re1 ated to the value of the incident angle. The software does 

not produce this calculation. For further discussion see Appendix A. 
The software has been applied since early June 1984 at NASA Langley f, 

Research Center. An average of ten data files have been stored daily, 

about half of which indicate plane wave propagation with printed delays 

and angles. This printed data is then compared with the weather reports, obtained from the meteorology center at NASA, on a daily basis. 

Finally, the data files stored on floppy disks will undergo a comprehensive spectral analysis, vainly for signature identification and power 

law verification. The software used for the spectral analysis is, or 

has stemmed from, DAISE (Digital Analysis of Infra-Sound Experiments) 

[23], that was developed previous to the course of this research. 

Daily evaluation of the data proved that most files stored by the 

system can be grouped into four categories, with their typical sources, 

as follows: 

*Within one mile from the center of the array. 

Man-made infrasounds were to be expected since the microphone array i s 

located in Langley Research Center where a wide range of activities take 

place on a regular basis. These include operating a1 1 sorts of machi nery, wind tunnels, et~. Future plans for near-field (both categories) 

signal analysis exist but they are not of interest at this point. 

Infrasounds produced by far-field events are discussed in the next two 

Weather-re1 ated 

IY~ crobursts 

Severe storms 

Near-f i el d* 

Man-made 

He avy machi nery 

Far-f i el d Mi 1 i tary aircraft s 
sections. The man-made category is discussed in section 4.4 while the 

weather-re1 ated category is discussed in section 4.5. 

4.4 Man-Made Si gnal s 

This section discusses the man-made signals with long-range propagation detected by the automatic detection routine described in the 

previous section. This category of infrasounds is of interest to this , 

study for two reasons. First, it verifies the capability of the system 

to detect infrasonic signals and to locate their sources. Second, it 

illustrates the necessity of more sophisticated signal processing techniques to distinguish between these events and tne meteorological ones. 

Among others, the non-weather re1 ated signal s detected include 

single-component infrasounds (for example 7.2 Hz, 10 Hz, 16.8 Hz), sonic 

booms generated by supersonic vehicles, and i nfrasounds from jet powered 

aircraft. Single component infrasounds have been detected daily. The 

system estimate of their directions places their sources in the main 

research center of Langley. .Thus, it is be1 i eved that they are mostly 

related to the different activities in the Center. Few sonic booms were 

detected during the summer of 1984. They were mainly associated with 

flights of supersonic military aircraft such as the Air Force F-15 which 

are stationed nearby at Langley Air Force Base. The high altitude of 

these sonic booms was indicated by elevation angles of 25" to 40" evaluated from the data. The third type of man-made signal was found in 

several files over a period of about three months. The common characteristics between tnese infrasounds include a sound pressure level of 
about 90 dB, frequency band of about 1 to 6.5 Hz, and a source direction 

ranging from 130" to 150°, re1 ative to true north. The consistency of 

these values between files, the relatively strong amp1 itude, and the 

directions found all have led to the belief that these infrasounds are 

generated during the take-off of mi 1 i tary jet powered aircraft from 

Langley Air Force Base. The location of the runway at the base re1 ative 

to the microphone array <is shown in Figure 4.25. 

One of the jet data files is considered here for illustration. The 

file was stored on August 28, 1984, by the automatic detection routine. 

The parameters of the software and the hardware were set to their nominal values as specified in Table 3.1. The signal on one channel before 

and during the event is shown in Figure 4.26 in the time domain and in 

Figure 4.27 in terms of its auto-power spectrum. The signal on microphone 1, 2, and 3 during a different time frame, of the same file, is 

shown in Figure 4.28. To demonstrate the convergence of the algorithm 

for this type of data, the filter response was recomputed for C1 x C2 in 

Figure 4.29, C2 x C3 in Figure 4.30, and C3 x C1 in Figure 4.31. This 

response was recomputed by feeding alternative data points of the time 

histories to the adaptive filter. This effectively reduced the sampl ing 

rate to 25.6 Hz (51.2 Hz originally). A 25.6 Hz sample rate is adequate 

in this case since the signal has a narrow frequency band as can be seen 

in Figure 4.27. Keeping this rate in mind, examining the filter 

response yields the following results 

C1 lags C2 by 0.742 sec. (-19 resolutions) 
aIOIMAL PAGE IS 9 5 

QUALsW 

Figure 4.25. Location of the runway in iangley Air Force Base relative 

to the microphone array. 
I 1 I I I I I I 

0 2.5 5 7.5 10 

Discrete time (in seconds) 

discrete time (in seconds) 

Figure 4.26. Time histories (in equivalent arbitrary units) of infrasounds generated by jet-powered aircraft, (a) before, and 

(bj during the event. 
12.8 Hz 

Di scre te frequency 

12.8 Hz 

Di screte frequency 

(b) 

Figure 4.27. Auto-power spectrums (in arbitrary units) of infrasounds 

generated by jet-powered aircraft, (a) before, and (b) 

during the event. 
I I I. 1 I I 

0 2.5 5 7.5 10 

Discrete time (in seconds) 

I I I I 1 1 I 1 

0 2.5 5 7.5 10 

Discrete time (in seconds) 

I I I I I I I I 

0 2.5 5 7.5 10 

Discrete time (in seconds) 

Figure 4.28. Time histories (in equivalent arbitrary units) of infrasounds generated by jet-powered aircraft, sarnpl ed on (a) 

channel 1, (b) channel 2, and (c) channel 3. 
Number of iterations 

(a > 

Weight number 

(b 

Figure 4.29. Filter response for jet data, C1 x C2, (a)) TDE vs. 

time, (b) Final weight vector. 
Number of iterations 

Weight number 

Figure 4.30. Filter response for jet data, C2 x C3, (a) TDE vs. time, 

(b) Final weight vector. 
Number of iterations 

(a) 

Weight number 

(b) 

Figure 4.31. Filter response for jet data, C3 x C1, (a) TDE vs. 

time, (b) Final weight vector. 
C2 leads C3 by 0.625 sec. (16 resolutions) 

C3 leads C1 by 0.117 sec. (3 resolutions) 

Based on these de1 ays, a source direction of about 145" can be found. 

The estimate of the bearing places the source roughly on the middle of 

the runway as seen in Figure 4.25. This result among others obtained 

from similar files presents a reasonable assurance that the system is 

indeed capable of detecting and locating the sources of infrasounds. 

4.5 Weather-Re1 ated Data 

As was mentioned in Chapter one, the far-field meteoroloyical 

events of interest in this research are clear-air turbulence and severe 

storms. There has not been clear evidence of clear-air turbulence 

detection on the system since it has been in operation. It is also 

believed that only a few storms were recorded. The uncertainty governing their detection -is high in all cases except for one severe storm. 

This storm was the hurricane "OIANA." It was relatively close (200 

miles) and stationary for at least one day. The data files stored by 

the system during that day, and their evaluation, present enougn necessary and reasonably sufficient conditions to identify the data with the 

hurricane. An overall evaluation of the system, including a discussion 

of the various factors that might have contributed to the lack of meteorological data, is discussed in the conclusion chapter. This section 

will concentrate on the hurricane and the issues surrounding its detection. 
A tropical storm developed into a hurricane over the coastal water 

of North Carolina during the period of September 10 to September 12, 

1984. It reached its full strength by September 12 wit0 recorded wind 

velocity up to 110 mph. The location of its "eye" on that day was relatively stationary about 200 miles away from NASAILangley Research Center 

at 1 atitude 33.8 and longitude 77.3. During the afternoon and evening 

of the same day, the system recorded eleven files, four of which pointed 

toward the site of the hurricane with a calculated sound pressure level 

of 85 dB. The first three files, were stored within approximately a 40 

minute time frame. This raises the possibility of having all three 

files associated with the same event. Examination of the spectral shape 

and contents of the signal on one channel through all four files supports this possibil ity. For reference, information in the signal sampled on one microphone is presented in Figure 4.32 as time history, in 

Figure 4.33 as auto-power spectrum, and in Figure 4.34 as log-amp1 i tude 

vs. 1 og-frequency aut o-power pl ot s .* Compared to other data f i 1 es recorded, these data present three, somehow unique, features. First, the 

signal has a wide frequency band. This can not be seen on the linear 

scale of Figure 4.33, because of the difference in the dynamic range of 

the components on both ends of the spectrum. But it is clear on the 

log-scale of Figure 4.34 where the band plotted is 240-16 Hz. Second, 

the roll-off of the spectrum is consistent with the prediction of 

-3.5 Meecham and Ford [21]. They predicted a rol 1-off proportional to u , 

*Software used to obtain these plots was developed at NASAILangley Research Center, by Dr. A1 1 an J, Zuckerwar and David Katzoff . 

C," 

L 

522 
apnq L [due-60 

1 

h h 

ra v Y V 

3 

a,?" 

L aJ 

CLL 

oms 

0 
while the ones computed in Figure 4.34 range between w 

-2.41 and .-3.05 

Third, the frequency band as well as the roll-off were consistent 

through all the files even though the fourth file was sampled hours 

after the first three. 

To illustrate the consistency between the microphones outputs, the 

signal on all three microphones during the same time period was plotted 

in Figure 4.35. Its auta-power spectrum is plotted in Figure 4.36 on a 

linear scale and in Figure 4.37 on a log-scale.* The adaptive filter 

algorithm was reapplied to the data in Figure 4.35 to demonstrate the 

type of convergence obtained for this signal. The response is shown in 

Figure 4.38, Figure 4.39, and Figure 4.40 for C1 x C2, C2 x C3, and C3 x 

C1, respectively. As can be seen in the figures, the TDE vs. time plots 

reflect the noisy (actually windy) environment in which the signal traveled and was then sampled. The directions evaluated from the different 

four data files were, with respect to the bearing resolution, within an 

acceptable tolerance from each other. A mean direction of 187", from 

north, was obtained. Tiis bearing is 'to be compared with the actual 

value of 193" as illustrated in Figure 4.41. 

The foregoing analysis and results clearly constitute the necessary 

conditions for true detection of the hurricane. The sufficiency, however, is lacking a higher detection rate, on behalf of the system, over 

a longer period of time. This issue, among others, will be di scussed in 

Chapter 5. 

*Software used to obtain these plots was developed at NASA/Langley 

Research Center, by Dr. Allan J. Zuckerwar and David Katzoff. 
'. . I:. ... -:.:.::: .:. 

3 ,.: ... 

.r 

Discrete time (in seconds) 

Figure 4.35. Time hi stories (in equivalent arbitrary units) of 

hurricane signal simultaneously sampled at (a) channel 

one, (b) channel two, and (c) channel three. 
Figure 4.36. Auto-power spectrums (in arbitrary units) of hurricane 

signal simultaneously sampled at (a) channel one, (bj channel two, and (c) channel three. 

\ 
a 

-0 

3 

C, 

SF- 

? 

P 

5 I 

U-J 

0 

-I 

Figure 4.37. Auto-power spectrums on a log scale (log-amp1 itude VS. 

log-frequency, in equivalent arbitrary units) of 

hurricane signal simultaneously sampled at (a) cnannel 

one, (b) channel two, and (c) channel three. 
Number of iterations 

(a> 

Weight number 

(b 

Figure 4.38. Filter response for hurricane data, C1 x C2, (a) TDE vs. 

tine, (b) Final weight vector. 
Number of iterations 

(a) 

Weight number 

(b) 

Figure 4.39. Filter response for hurricane data, C2 x C3, (a) TDE vs. time, 

(b) Final weight vector. 
Number of iterations 

(a) 

Weight number 

Figure 4.40. Filter response for hurricane data, C3 x C1, (a) TDE vs. 

time, (b) Final weight vector. 
mlGBNAL PAGE bS 

OF 114 

Figure 4.41. Hurricane site relative to the microphone array in 

NASA Langley Research Center. 
4.6 Summary 

0 

The response of the basic LIYSTUE algorithm was evaluated for deterministic data sequences. The simulation justified the specific choice 

of the filter parameters. The effect of varying the signal characteristics on the response was demonstrated. A direct proportionality between the SNR and the rate of convergence of the weight vector was observed. The dependency qf the main lobe width of the weight vector on 

the frequency band of the signal, as predicted by Ahmed and Carter [17], 

was verified. The ability of the algorithm to track time-varying delays 

was also illustrated in the simulation. A steady state time-lag between 

the actual varying delay and the estimate of about 64 samples was observed for SNR=2. The algorithm was then tested using Space Shuttle 

VIII launch noise. The response was satisfactory in the sense that tne 

1 aunch site was pi n-poi nted within the bearing resolution. The specific 

software used in the field on a daily basis was described. It consisted 

of an automatic detection and bearing eval uati on routine. The operati on 

of the system over the period of June to September 1984 was evaluated 

with respect to the rate and type of signals detected. Data evaluation 

proved that man-made signals dominated the data files stored by the 

system. A sample of detected signals from jet-powered aircraft was 

presented. This sample verified the capability of the system to detect 

infrasounds and locate their sources. The system detection of weatherrelated infrasounds was not conclusive. Only one severe storm was detected with reasonable assurance. The assurance of the detection was 

based on the consistency of the signal characteristics and bearings 

through several files sampled on the day of the storm. 
CHAPTER 5 

CONCLUSION 

5.1 Remarks 

The goal of this research was the design of a system for detection 

and bearing eval uati on of i nfrasounds. Based on the resul ts presented 

in Chapter 4, it is be1 ieved that this goal has been achieved. The 

adaptive LMSTDE algorithm employed in the design was evaluated in section 4.1. Simulation results were in total agreement with theoretical 

expectations. The response of the algorithm to the shuttle noise clearly demonstrated its ability to estimate delays in cases of long-range 

infrasonic propogation. The automated-system performance, in the actual 

field, proved effective through detection and bearing evaluation of manmade signal s. This class of signal s consisted of several types of 

infrasounds incl uding that of jet powered aircraft, whose signature i s 

similar to what is anticipated in meteorological events. Acceptable 

performance of the system, as an infrasonic detector and bearing evaluator, was a1 so achieved for 1 arge-scale weather events such as the 

hurricane "DIANA." 

This work, however, is faced with the difficulty of evaluating the 

system as a turbulence detector for two reasons, both of which are 

beyond the designer's control. First, the data obtained in this case is 

poorly defined. That is, it is not repeatable, nor is it unambiguously 

verifiable. The second reason is the 1 ack of a 1 arge nunber of weather- 
related data files needed for the evaluation. As pointed out earlier, 

few meteorological events were detected on the system. Only the hurricane was identified with reasonable assurance. Some of the factors that 

mignt have contributed to the low detection rate of events, or of files 

within the same event, are listed below. 

1. Low sound pressure level of weather-re1 ated infrasounds. A 

minimum level of 54 dB, derived from the system calibration, is 

required at the site of the array for the signal to be detected. The fact that the hurricane, despite being a powerful 

event, produced only 85 dB, raises the question of how big and 

how close the event should be in order to be detected. This 

question is left open for acoustical investigators. 

2. Noisy environment. The presence of strong man-made signals at 

the array location makes it difficult for the system to identify smal 1-scale eath her events. 

3. dind effects. Small wind fluctuations at tne array can be 

considered uncorrelated noise and usually do not affect tne 

system response. However, this is not the case in windy conditions, Strong wind may correlate at two or more microphones, 

destroying the coherence property of the turbulence signal. 

4. Mu1 tiple del ay. In situations of scattered thunderstorms or 

large frontal systems, the signal sampled at the array is the 

superposition of infrasounds from more than one event. The 

multiple delay pattern embedded in this type of signal creates 

an ambiguity problem that can not be solved on the existing 
system. <! 

5. Short-term stationary signals. A1 though the lifetime of a 

meteor01 ogical event i s re1 at ively 1 ong (several hours i n some 

cases), the period over which the turbulence, and hence the 

infrasonic emission, exists is virtually unknown. Therefore, 

the event may be easily missed in the case of turbulence that 

is stationary ov,er a short time relative to the computational 

speed of the system; i .e. 20 seconds worth of data are sampled 

every 7 minutes (average computational time), so infrasoni c 

emission for less than 7 minutes can be missed. 

5.2 Future dork 

More work is needed to improve the system performance for turbulence identification purposes. One work area is the investigation of 

signal processing techniques that are appropriate for el imination of 

non-weather data and removal of wind effects from the signal before 

applying the adaptive LMSTDE a1 gori t hm. Once enough know1 edge i s 

obtained about the infrasonic signature of turubulence, pattern recognition methods can be used to discriminate against non-weather data. 

Removal of wind effects may require the development of a wind-measurement system that is synchronized with the infrasonic detection scheme. 

The purpose of this system would be to provide the wind information 

needed for compensation in the system input. 

Another future research consideration is. the study and modification 

of the LMSTOE algorithm for supressing the side lobes in the weight 

vector and processing of multiple delay signals. Both issues are 
important for avoiding any ambiguity in the weight vector peak that may 

be caused by the characteristics of the input signal. 

Future work also needs to consider implementing the LMSTDE algorithm in hardware using dedicated signal processing chips. Hardware 

implementation would dramatical ly improve the computational speed, and 

hence, the detection rate, of the system in tne case of short-term 

stationary signals. t 
REFERENCES 

"Low-Altitude Wind Shear and Its Hazard to Aviation," by National 

Research Counci 1 ., National Academy Press, Washington, D.C., 

1983. 

Airline Pilot (Technical Talk), Vol. 46, No. 4, April 1977. 

P. Kuhn, Clear Air Turbulence: Detection by Infrared Observations of Water Vapor Science, 196, pp. 1099-1100, 1977. 

Allan J. Zuckerwar and Harlan K. Holmes, "A Unified Acquisition 

System for Acoustic Data." NASA Technical Note TN 0-8327, March 

197 7. 

C. ti. Knapp and G. C. Carter, "The Generalized Correlation Method 

for Estimation of Time Delay," IEEE Trans. Ascoust., Speech, 

Signal Processing, Vol. ASSP-24, No. 4, August 1976. 

J. C. Hassab and R. E. Boucher, "Upti~nurn Estimation of 'Time delay 

t3y A General i zed Correl ator ," IEE Trans. Acoust . , Speech, Si gnal 

Processing Vol. ASSP-27, No. 4, August 1979. 

A. Papoul i s, Probabi 1 ity, Random Variables and Stochasti c 

Processes. New York: McGraw-Hi 11, 1965. 

G. M. Jenkins and 0. G. Watts, "Spectral Analysis and Its 

Appl ications ," San Francisco, CA, Holden-Day, 1968. 

P. R. Roth, "Effective Measurements Usin Digital Signal Analysis," IEEE Spectrum, Vol. 8, April 1971. 

G. C. Carter, A. H. Nuttall, and P. G. Cable, "The Smoothed 

Coherence Transform (SCOT) ," Naval Underwater Systems Center, New 

London Lab, New London, CT, Tech. Memo TC-159-72, August 8, 1972. 

B. Widrow, "Adaptive Filter, in Aspects of Network and System 

Theory," R. Kalman and N. DeCl aris, Eds., New York: Holt, 

Rinehart and Winston, 1971, pp. 563-587. 

Steven A. Tretter, Introduction to Ui screte-Time Signal 

Processing," New York: John Wiley and Sons, 1976. 
L. J. Griffiths, "Rapid Measurement of a Digital Instantaneous 

Frequency,I1 IEEE Trans. Ascoust. Speech, Signal Processing; Vol . ASSP-23, Apri 1 1976. 

B. Widrow, et a1 . , "Adaptive Noise Cancelling Princi ples and 

Appl ications ,I1 Proc. IEEE Vol . 63, December 1975. 

8. Widrow, P. Montey, L. Griffiths, and B. Goode, "Adaptive Antenna Systems," Proc. IEEE, Vol. 55 

F. A. Reed, P. L. Feintuch, and N. J. Bershad, "Time &lay tstimation Using the LMS Adaptive Filter Static Behavior," IEEE 

trans. Acoust. Speecn, Signal Processi ng, Vol . ASJP-29, June 

1981. t 

O. H. Youn, N. Ahmed, and G. C. Carter, "An Adaptive Approach for 

Time kl ay Estimation of Band Limited Si gnal s,I1 IEEE Trans. 

Acoust., Speech, Signal Processing, Vol. ASSP 31. 

P. L. Feintuch, N. J. Bershad, and F. A. Reed, "Time Delay Estimation Using the LMS Adaptive Fi lter-Dynamic Gehavior," IEEE 

Trans. Acoust. Speech Signal Processing, Vol . ASSP-29, June 1981. 

John A. McDonal d, "Natural ly Occurring Atmospheric Acoustical 

Signals," J. Acoust. Soc., Vol. 56, No. 2, August 1974. 

Eric S. Posmetier, "1-to-16-Hz Infrasound Associ ated with Clear 

Ai r Turbulence Predictors," Journal of Geophysical Research, Vol . 79, No. 12, April 20, 1974. 

W. C. Meecham and G. W. Ford, "Acoustic Radiation from Isotropic Turbulence," The Journal of the Acoustical Society of Pmerica, 

Vol. 30, No. 4, April 1958. 

ivl. Tahira and N. L. Donn, "Anomalous Infrasound from Space 

Shuttle and Skylab I," J. Acoust. Soc. h., Vol. 73, No. 2, February 1983. 

C, 5. Khalaf and J. k. Stoughton, "Software tlevelopment for 

Infrasound ivleasurernent System," prepared for NASA Langl ey 

Research Center, under NAS1-17099-Task No. 10, August 1983. 
APPENDIX A 

Thi s appendix represents the derivation of a re1 ationshi p between 

the source bearing and source elevation, and the delays obtained for an 

arbitrary sound wave. A general relation is first derived for one pair 

of microphones. The specific re1 ations for the array used in this work , 

are then discussed. 

Figure A.l illustrates the elevation effects on the delay between 

two microphones, where 

8 is the source direction relative to the line joining the 

microphones 

4 is the source elevation angle 

D' is the delay due to the el evated-source signal 

D is the projection of D' on the horizontal plane (delay for $=0) 

V is the sound velocity 

2 is the distance between microphones. 

Inspection of the figure yields the following two relations, 

D'xV - Dl cos 4 = - - - 

DxV D 

and Dx V 

cos 8 = - 
Figure A. 1. Elevation effect. 
substituting for D in (A.l) gives 

V D ' COS$ =- xR cos 8 

Equation (A.3) is a transcendental equation that re1 ates the direction 

and elevation of the source to the delay between two microphones. Note 

that Dt is given, as the output of the LMSTDE algorithm, while 4 and 

e are to be found. 

The array configuration, and the values of 8 (el, e2, 6,) for the 

three pair of microphones, are shown in Figure A.2 (horizontal plane). 

I I I 

Let 4 , C$ , and I+ be the del ays obtained from micropilones (3, 4), 

(3,2), and (3,1), respectively, due to a signal propogating with 4" 

elevation. Applying equation (A.3) to the different pairs yields 

v 

cos $ = - D2 

R cos 82 

I 

I 

v 

cos 4 = - D3 

R COS 83 

(A. 4c) 

el, 82, and 83 represent the source direction relative to the 

microphones axis as shown in the figure. If B represents the source 

direction re1 ative to true North, then 81, 82, and 83 can be expressed 
North I 

Figure A.2. Relative angles. 
in terms of $ by 1 inear relationships, i .e., in figure A.2, the re1 ations are 

, I I I 

Hence, having 4, 4, and 4 estimated by the system, the problem of 

finding 4 and B becomes that of solving two equations (out of the 

three in eq. (~.4)) with two unknowns. 

Note that the re1 ations of (A.5) hold for a certain range of B 

only. Thus, when solving for the elevation and the exact bearing of the 

source in a given situation, these re1 ations have to be evaluated with 

the bearing values provided by the system.