Group Title: signal detection from noise utilizing zero-crossing information
Title: A signal detection from noise utilizing zero-crossing information
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00098936/00001
 Material Information
Title: A signal detection from noise utilizing zero-crossing information
Physical Description: vi, 136 leaves : ill. ; 28 cm.
Language: English
Creator: Yeo, Woon Cheon, 1944-
Copyright Date: 1975
 Subjects
Subject: Signal theory (Telecommunication)   ( lcsh )
Random noise theory   ( lcsh )
Matrices   ( lcsh )
Electrical Engineering thesis Ph. D
Dissertations, Academic -- Electrical Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
 Notes
Statement of Responsibility: by Woon Cheon Yeo.
Thesis: Thesis--University of Florida.
Bibliography: Bibliography: leaves 132-135.
General Note: Typescript.
General Note: Vita.
 Record Information
Bibliographic ID: UF00098936
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000171239
oclc - 02954412
notis - AAT7664

Downloads

This item has the following downloads:

signaldetectionf00yeow ( PDF )


Full Text










A SIGNAL DETECTION FROM NOISE
UTILIZING ZERO-CROSSING INFORMATION






By

WOON CHEON YEO


A DISSERTATION PRESENTED TO TEF GRADUATE COUNCIL OF
THE UNI'VEhSITY CF FLORIDA IN PARTIAL
FULFILLUENTO OF THIE REQUIREMENTS FOR THE DEGREE OF'
DOCTOR CF PHILOSOPHY




UNIVERSITY OF FLORIDA
1975


























UNIVERSITY OF FLORIDA
11111111111111111111111111111111111111
3 1262 08666 219 3 1





























Copyright by
Woon Cheon Yeo
1976














ACKNOWLEDGMENT


Regretful to leave the University of Florida

Where knowledge is delivered with such warmth 'n

Where friendly hearts always surround you.



Supported and given the opportunity to delve into the

Meaning of knowledge,

Grateful to Dr. Jack R. Smith as

I owe all these to his generosity.



Admire and take pride in all my committee as

They embed confidence in your strength.














.TABLE OF CONTENTS


Page

ACKNOWLEDGMENT . . . . . . . . .. iii

ABSTRACT . . . . . . . . ... .. ... . v

CHAPTER I. INTRODUCTION . . . . . . ... 1

CHAPTER II. ZERO-CROSSING INTERVAL DENSITY FUNCTIONS 8

Derivation of Zero-crossing Interval Density
Function . . . . . . . . . . 8
Evaluation of W for Gaussian Process . . . 12
Zero-crossing Interval Density Function for a
Bandlimited White Gaussian Noise . . . . 20
Evaluation of W for Sine Wave Plus Noise . . 23
Zero-crossing Interval Density Function for
Sine Wave Plus Noise . .. . . . . 30

CHAPTER III. IMPLEMENTATION AND EVALUATION OF THE
DETECTOR . . . . . . . . ... .. 38

Implementation of the Detector . . . . . 38
Comparison to Correlation Detector . . . 54

CHAPTER IV. FURTHER DISCUSSIONS AND CONCLUSION . . 65

Operating Characteristics in Practical Situations 65
Conclusion . . . . . . . .. 75

APPENDIX A. JACOBI'S THEOREM . . . . .... 79

APPENDIX B. EVALUATION OF A DETERMINANT . . ... 83

APPENDIX C. NUMERICAL INTEGRATION . .. . ... 87

APPENDIX D. INTERPOLATION BY SPLINE FUNCTION ... . 90

APPENDIX E. LIST OF PROGRAMS . . . . . . 93

BIBLIOGRAPHY . . . . . . . . . . 132

BIOGRAPHICAL SKETCH . . ... . . . 136














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



A SIGNAL DETECTION FROM NOISE
UTILIZING ZERO-CROSSING INFORMATION


By

Woon Cheon Yeo

December, 1975


Chairman: Dr. Jack R. Smith
Major Department: Electrical Engineering


A method of detecting signals from noise utilizing the

zero-crossing information is analyzed.

When the noise is stationary white Gaussian, a corre-

lation detector is considered to be optimal. But, if the

process of signal and noise is complex, then a nonlinear

detector may give a superior result in overall performance.

In a situation where the signals and noise have

nonstationary statistics, a detector which measures the

zero-crossing intervals of the process has been found to be

very efficient in detecting signals from noise. The operating

characteristic of the zero-crossing measuring detector is

analyzed in a model with signals of sine wave form imbedded

in a stationary bandlimited white Gaussian noise. Then, to








observe the operating characteristic of the detector in a

complex situation, the detector is analyzed under deviated

conditions of signals and noise from the assumptions. The

performance of a correlation detector is also analyzed under

identical conditions, as a comparative reference to the

performance of the zero-crossing measuring detector.

The operating characteristic of the zero-crossing

measuring detector is shown to be effective in detecting

phasic events from a sleep electroencephalogram, where the

process of the signals and noise is very complex.














CHAPTER I
INTRODUCTION



An optimum detector of a signal from noise may well be

described as "a detector which best satisfies given criteria

under a given set of assumptions" (Whalen).

The first problem in designing an optimum detector is,

then, making proper assumptions on parameters that completely

describe the process of signal and noise, and making proper

selections of parameters that are relevant to discriminating

signals from noise. If either the real process is different

from that described by the assumptions on parameters, or any

of the relevant parametersare neglected, then the detector

may be only an approximation to the optimum.

In the case when the process is complex, we may not be

able to make assumptions on parameters, or even if we could

define the process by the parameters, an optimum detector

based on all the relevant parameters becomes too complicated

to implement. Thus we need to select among the definable

parameters the most efficient features from the process as

the discriminating parameters so that a detector based on

these parameters should give an acceptable performance in

the complex situation. The efficiency of the parameters is

judged by the given criteria under the assumptions on the

parameters. The performance of the detector based on these








efficient parameters, however, may not be optimum if the

environment of signal and noise changes so as to include

other efficient definable parameters.

In a stationary Gaussian noise, a correlation detector

is considered to be optimal (Wainstein). If the noise process

is complex so that it could not be simply definable as sta-

tionary Gaussian, then a superior detector may require a

nonlinear characteristic, and the implementation of the

detector may be complicated and costly. Especially if the

process is not stationary, then the adaptation of the detec-

tor becomes a serious problem, and it may leave the detector

complicated to operate. Thus, in general, the design of a

detector should include considerations such as performance,

cost of implementation, simplicity of implementation, and

simplicity of operation. The priority of the considerations

may vary depending on the circumstances of the particular

application of the detector.

In designing a system to detect phasic events from

sleep electroencephalograms(EEGs), one encounters the

following complexities in the process of the EEG activities:

The phasic events, which are designated by arrows in Figure

23, are distinguished from the irregular background EEG by

their resemblance to sinusoidal waves. Sigma spindles, for

example, last from 6 to 20 cycles with a frequency between

12 and 15 Hz depending on the individual. The amplitude of

spindles and background EEG activities, which is considered

as noise when we desire to detect sigma spindles, also








varies widely depending on the individual and the electrode

resistance. The a priori probability densities of the fre-

quency and amplitude distribution are not usually known.

Even though the strength of the signal (sigma spindles) and

noise varies, the signal-to-noise ratio remains relatively

constant; that is, a subject with high background EEG acti-

vity also has high amplitude in the spindle activity. The

strength and the power spectra of the background EEG undergo

a continuous change as a sleep progresses.

In practical applications of the phasic event detector,

it is desirable that the detector be able to perform equally

well an any subject without a priori information on the

frequency and the strength of the EEG activity. To accomplish

this a detector could be designed which adapts to the varying

frequency and intensity of the EEG activity, but it might be

too costly to implement.

Smith et al.(43,44) designed a phasic event detector

which employs the zero-crossing interval of the EEG as the

discriminating parameter, and they obtained rather superior

performance from the detector. The detection rate was high

enough to meet the requirement, while the false alarm rate

was extremely low for any kind of noise situation. The

detector could be applicable to a wide variety of subjects

without any adjustment of thresholds to accommodate each

subject. It was simple to implement, and its overall perfor-

mance was superior to a correlation detector based on weak

assumptions on the signal and noise. And thus a mathematical


__








analysis of the zero-crossing measuring detector was

suggested and carried out in this study to provide a better

understanding on the efficiency of the feature of zero-

crossing interval as a discriminating parameter of signals

from noise.

Since both the signal(phasic events) and the noise(back-

ground EEG activity) processes in sleep EEG are too compli-

cated to be explicitly defined, an analysis of any phasic

event detector from sleep EEG is necessarily based on improper

assumptions. Also any attempt to make the model close to the

reality would leave the analysis untenable. Thus the analysis

of the detector is carried on with an idealized model of

signal and noise. The signal is modeled by sinusoidal waves

and the noise is assumed to be bandlimited stationary Gaussian.

Then, to observe the performance of the detector in a nonideal

situation, the operating characteristic of the detector is

analyzed under conditions of signal and noise different from

the assumptions. In each case the performance of the correl-

ation detector is also analyzed as a comparative reference to

the zero-crossing measuring detector,

To compute the detection and false alarm probabilities of

the zero-crossing measuring detector, we need the zero-cross-

ing interval density functions of the processes with noise

alone and signal plus noise. The distribution functions for

the successive zero-crossing intervals of the bandlimited

Gaussian noise have been developed by Rice(36,37), McFadden

(26,27) and Longuet-Higgins(23,24). Experimental works on


S.








the zero-crossing interval distribution of Gaussian noise

have been presented by Rainal(32) and Blotekjaer(6).

Rice's zero-crossing interval density function agrees well

with experimental data for smallvalues of zero-crossing

interval. The zero-crossing interval density function derived

by McFadden gives good approximations to Gaussian processes

having arbitrary power spectra. For the particular Gaussian

process which has an ideal low pass spectrum, the approxi-

mation by Longuet-Higgins fits closely to the experimental

data given by Blotekjaer. Longuet-Higgins derived a proba-

bility density function of the spacing between the ith zero

and the (i+m+l)th zero of a stationary random function. The

density function is expressed as a rapidly convergent series,

and, when the process of the random function is Gaussian, the

first two terms of the series may be expressed in terms of

known functions. This Longuet-Higgins approximation is used

here to compute the zero-crossing interval distribution of a

bandlimited Gaussian process.

In the signal-plus-noise case, the zero-crossing inter-

val density function is approximated by the first term of

the Longuet-Higgins series. This approximation is also

derived by Rice(38) for Gaussian processes. Cobb(lO) gave a

treatment on this approximation for a process of sine wave

plus bandlimited Gaussian noise. He reports in the same paper

about the zero-crossing interval measuring principle in

distinguishing two signals separated in frequency. Yet his

paper does not provide details of the detector and any








comparative result to other kinds of detecting methods. The

computation of the zero-crossing interval density for sine

wave plus noise in this study is based on the treatment of

Cobb.

Derivations of the zero-crossing interval density func-

tion for Gaussian noise and Gaussian noise plus signal are

presented in Chapter II. An expression of the first two terms

of the Longuet-Higgins' series in terms of the autocorrelation

function of the Gaussian noise is derived. The computed

density function is plotted along with the experimental data

as a comparison. The derivation and the computed data of the

zero-crossing interval density for signal plus noise are also

presented in Chapter II.

In Chapter III, the schematic of the zero-crossing

measuring detector is described and the method of maximizing

the detection probability for given false alarm rate is

discussed. The performance of the detector is compared to

that of the conventional correlation detector under the iden-

tical situation of signal and noise.

The advantages and disadvantages of the zero-crossing

measuring detector compared to the correlation detector are

further discussed in Chapter IV, with practical considerations

in detecting phasic events from sleep EEGs.

All the programs to compute various functions were

written in FORTRAN and assembly language, and they were

processed through a minicomputer PDP-8/e. Without an extended

arithmetic unit, the speed of the PDP-8/e was very slow to

compute some numerical integration. Thus, to save the


1_ _






7

computing time, we employed data interpolations in the

computation of the zero-crossing interval density functions

and evaluated the statistical functions from their series

expansions.













CHAPTER II
ZERO-CROSSING INTERVAL DENSITY FUNCTIONS


Derivation of Zero-crossing Interval Density Function


For simplicity, we will state that a random function
f(t) has a zero or a crossing at t=t1 if f(t1)=O. Similarly

we use up-crossing or down-crossing to indicate a crossing
with f'(t,)>0 or f'(t1)<0 respectively.

Consider the instance when a random function f(t) has
up-crossings in the infinitesimal intervals (ti,ti+dti)
(i=l,...,n), and let us denote the joint probability of this
happening by W(+,+,...,+)dtldt2...dtn. To denote the joint
probability that the random function has down-crossings in
some intervals, we change the plus sign to minus at the
corresponding positions in W. Thus W(+,-,+)dtldt2dt3, for

instance, would represent the joint probability that a random
function f(t) has an up-crossing in the interval (tl,t.+dtl),
down-crossing in (t2,t2+dt2), and an up-crossing in

(t3,t3+dt3).
We next introduce p(m;r) to denote the probability
density of the interval between the ith and the (i+m+l)th
zeros of f(t). This probability density is denoted by Pm(r)
by other authors, but, aince in this paper a probability is

represented by an uppercase P and a probability density is








written in a lowercase p, the symbol p(m;r) is used here.
Then p(O;T)dr is the probability that the interval between

two successive zeros has length r, and p(Ojr)drW(+)dt1
expresses the probability that f(t) has an up-crossing at
t=t1 and a down-crossing after an interval r without having

any other crossings in between the two crossings. Likewise
p(2n;r)drW(+)dt1 represents the probability of f(t) having
an up-crossing at t=tl and a down-crossing after an interval
r, containing 2n crossings anywhere in between the two
crossings.
If f(t) has an up-crossing at t=t1 and a down-crossing
at t=t2, then f(t) may contain 2n (n=0,1,2,...) crossings
between t=t1 and t=t2. W(+,-)dtldt2, the probability of f(t)
having an up-crossing at t=t1 and a down-crossing at t=t2,
is then the sum of the probabilities of all possible cases
when f(t) has an up-crossing at t=t1 and a down-crossing
after an interval T=t2-1t, including 2n (n=0,1,2,...) cross-

ings within the interval r. That is

W(+,-)dtldt2=[p(O;r)+p(2;-)+p(4;r)+...]drW(+)dt1, (1.1)

or, since dT=dt2,

W(+,-)/W(+)=p(0;r)+p(2;T)+p(4;r)+..., (1.2)

(McFadden 1958).

In a similar fashion we can derive an expression

related to W(+,-,-). Let us consider the situation when f(t)

has an up-crossing at t=t1 and a down-crossing at t=t ,







containing a down-crossing at t=t2 anywhere in between t1 and
t Since the crossings at t=tI and t=t3 are an up-crossing
and a down-crossing respectively, f(t) can only have an even
number of crossings in the interval (tl,t3), and one of the
down-crossings among these crossings should pass zero at t=t2.
Suppose we have 2n crossings in (tl,t3) where t3-tl=tr Then
the down-crossing at t=t2 could be one of those n down-
crossings, and, therefore, there are n different possibilities
of f(t) having a down-crossing at t=t2. Since the down-cross-
ing at t=t2 could be anywhere in (tl,t3), if we integrate the
probability of each possibility with respect to t2 from tl to
t31we get p(2n;r)dTW(+)dt1, and we obtain for the n possibili-
ties a total probability of np(2n;-)dTW(+)dt1. And the sum of
the probabilities for all possible values of n,

[p(2;r)+2p(4;r)+...+np(2n;r)+... drW(+)dtl, (1.3)

is then the probability of f(t) having an up-crossing at t=t1
and a down-crossing at t=tl+T containing at least one down-
crossing at t=t2 anywhere in between tI and t3.
Now if we integrate W(+,-,-)dtidt2dt3 with respect to t2
over the interval (tlt3), we obtain the probability of f(t)
having an up-crossing at t=tI and a down-crossing at t=t3
containing at least one down-crossing at t=t2 anywhere in
between tI and t3. This probability should be equal to that
expressed by Equation (1.3). Therefore we obtain (Longuet-
Higgins),

S dt2 [W(+,-,-)A(+)] =
ti t2






By rearranging Equations (1.2) and (1.4),

p(o;r)=W(+,-)/W(+)-[p(2;T)+p(4;T)+...], (1.5)

p(2;T)= S dt2W(+,-,-)/W(+)] -
tl [2p(4;T)+3p(6;r)+...]. (1.6)

Substituting Equation (1.6) into Equation (1.5),

p(0;-)=W(+,-)/W(+) S dt2[W(+,-,-)/J(+)] +
t1 p(4;-r)+2p(6; )+3p(8;r)+.... (1.7)

If we neglect the terms [p(4;r)+2p(6;T)+3p(8;r)+...] in the
above equation, we get an approximate expression for the
zero-crossing interval density (Longuet-Higgins):

p(0;r)=W(+,-)/W(+)- f dt2[W(+,-,-)/W(-)]. (1.8)
tl

By neglecting the remaining terms in Equation (1.7) we
are ignoring the probabilities of f(t) having 4 or more
zeros within the interval of length T.








Evaluation of W for Gaussian Process


Let us consider the probability W(+,+,...,+)dt dt2...dtn

that f(t) has an up-crossing with an arbitrary positive slope
f'(t) in each of the small intervals (ti,t.+dti) (i=1,...,n).
For convenience write

f(ti)=fi,

f'(ti)=g

and let p(fl,...,fng'1.. ,gn) denote the joint probability
density of the fi and gi (i=l,...,n). Thus p(fl....,f ,
gl,...",n)df ...dfndgl ...dgn is the probability that the fi
and gi lie in the intervals (fi,fi+dfi), (gi,gi+dgi) respec-
tively for each i.
If f(t) has a zero in (ti,ti+dti) with a gradient gi,
then f(ti) must lie in a small range of value of extent

gil dti. Especially if f(t) has an up-crossing in the
interval (ti,ti+dti), then the range of the gradient is
O of p(fll ... f l ... ,gn) over the ranges (O
(-gidti up-crossings at the intervals (ti,ti+dti) (i=1,...,n). That
is,
W(+,...,+)dtl...dtn
00 00 0 0
= Odgl...dg S dfl... r df p(fl,..f, )
0 -(dt gtn1)
(2.1)







Since the range of f. (-g.dti.,) is small for each i,

we may consider P(fli...fng 1, ... g n) to be constant with

respect to f. within each interval (ti,t.+dt.). Then the

integrals in Equation (2.1) could be approximated as

W(+,...,+)dtl...dtn
CO 00 0 0
=Sdg ...SdgnP(f ... ,fn,g, .... gn ) d df ,
0 0 -gldt1 -gndt

which becomes

W(+,...,+)dtl...dtn

=fdgl... dgngl...gnp(...,0,gl,...,gn)dtl...dtn,
O 0

and thus,

W(+,.. ,+)
00 00
=.dgl... dgngl...gnp(0....,0,gl,...,g ). (2.2)
0 0

Now we try to find the expression of p(fl,...,f n

gl..." n) in terms of the autocorrelation function of the
Gaussian process f(t). The covariance matrix of the 2n

variables fl, ...,fngl...,gn is

M= Eflfi,... ,Eflfn,Efg ,...,Eflgn


EfEf EfEfngn (2.3)
Eglf,...,Eglfn,Egg ... ,Eg gn


S Egn fle...,Egn n,Egng],...,Egngn ,

where Efigj represents E[figjj, an expectation of figj.





14

Let R.. denote the autocorrelation function of f(t)
.J10
such as
Rij=R..ji=R(t-ti)=R(-)=E[ffif], (2.4)

where r=t.-t.. Then,
J 1

E[fi j]=E[f(t)-f(t)i ]
-j i t (tlt-'

1
=E[f (i)f(t+r) t=t

=E
= E[f(ti)f(ti+-)]

=6-R(T)=Ri (r), (2.5)

where the prime denotes differentiation with respect to f,
and 6 denotes partial differentiation.

E[gif ]=E[if(t) f(t )]

=E[ -f(t--r) f(t-)l

=Ev[6-f (tj-T)f(tj)]

6E[f(t -T)f(t.)]


=-R(-7)=-R(t)=- '(), (2.6)


EEgi g]=E[- f(t)| f(t)(;





1i


=ELif (t)l tti f(t+ T) I t=.t
dt It"^=t b

=E[d (t)=t f(t +r)]
6-r i
=-E[ F-Af(t)ttif(ti +1]
S'-dt t=t

= -EjE[ f (ti )f (t.i+) ]


=-6 R()=-R "(7).' (2.7)

Thus the covariance matrix in terms of R(r) is

M= R11,....,Rin, R11',"' Rn'
*-

Rnl ....,Rnn, I R R n'




-Rnl '***-Rnn''-Rnl ' nn

By the Gaussian hypothesis we have




=(2Tr)-nD exp(-iXTM-lX)
22n
=(2Tr)-rD fexp(- i ), (2.9)
i,j=l J1i;J )'

where X and XT are column and row vectors of 2n variables
(f', ... ,f,gl ..'.,gn) respectively, D=IMI, fi=gi' and qij
is the element of the ith row and jth column of the matrix
Q=M-1. We will denote a matrix with elements qij by (qij).





16


Substituting Equation (2.9) into Equation (2.2) we get

W(+,...,+)
i. 0 n
=(2w)-nD-~"dg ... Sdgn g2. .gnexp(- iE q nin+jgig ).
0 0 10=1
(2.10)
The summation in Equation (2.10) involves only the
last n rows and columns of (q..). We let (hij) denote the
inverse matrix of the matrix composed of the last n rows
and columns of (qij) such as


(hij)= qn+i,n+l"''qn+1,2n' -1
I I (2.11)
q2n,n+l" "' .2n,2n ,

and modify Equation (2.10) using this matrix.

W(+,...,+)



n
(2)--D-*(hij)l Jdgl.. dgn g2. gn Z(2,) ggn

(hj)l ~exp(-i, =1 n+in+j ig )


=(2Tr)- inDf i(h )I Sdgl... dg g g2...g Z(g,h), (2.12)

where

Z(gh)=(2w)-n(h)l-exp(-i 1 .n+in+jij)" (2.13.)

Now Z(g,h) is an ordinary normal distribution of n variables

(gl,..,g n) with the covariance matrix (hij). For conve-
nience, we normalize the covariance matrix (h i) such that
the (i,j)th element of the new covariance matrix is







1
vij=h i/(h iih. ),

and (v. ) is the covariance matrix of the new variables
wi=gi/hii2. Then, with the above change of variables,

Equation (2.12) becomes (Longuet-Higgins)

W(+,...,+)=(2T)-nD-1 (hij)1 (hlh22...hnn)aJ, (2.14)

where

J = Xdwl... dwnw w2"... nZ(w,v), (2.15)
0 0

and Z(w,v) is a normal probability density of n variables

(w.,...,w ) with the normalized covariance matrix (v ij).

By Jacobi's theorem (Appendix A), the determinant and

(i,j)th element of the matrix (hij) are computed as

I(h ij)I =D/D1, (2.16)

hij=D" R11...., R1n, R1j

(2.i7)
Rn...., Rnn, R
1' nn' nj
-Ril' ...'-Rin'-Ri"

where

Dl= R11,...,Rin

S (2.18)

R i, .... R .


Then, by substitution of Equation (2.16) into Equation
(2.14),
W(+,....,+)=(2)-InD -(hih22...hnn n. (2.19)





18

We now compute W for the special cases when n is 1, 2
or 3. If n=l, then Z(w,v).is a normal distribution of a
single variate, and we have
S001
Jl= Jdw(2T)-nwexp(-iw )=(2r)2,
0
D1=R11 =R(0),
and

h 1 R, R11' =R(O)-1 R(0) R(O)'
-R11',-R1" -R(0)',-R(0)"

=-R(0)".

Thus Equation (2.19) for the case n=l becomes

W(+)=(2T)-PD1-2h11J 1

=(2T)-1[-R(0)"/R(0)]P. (2.20)

When n=2 or 3, we may use the results of Nabeya(30)
and Kamat(19) to compute J2 and J3:

J2=(2y)- [(1-v12) 1Vcos1O(-12)], (2.21)

J3=(2n)-3/2[ (vij)lI+(slbi+s2b2+s3 3)], (2.22)

where
sl=cos- v31 v12-v23)(1-v312)-(1-v2 12

bl=v31v12+v23'

The angles of arc cosine are to be taken in the range (0,r),
and s2, s3, b2 and b3 are obtained by cyclic permutation of





19

the v... This gives

W(+,+)=(2T)-2 (R11R22-R12R21) (h11h22) (1-v122+

v12cos-1 (-12)], (2.23)
and
W(+,+,+)

=(2,1)-3-(hilh22h3 )[ (vj)i +(s1lbl+s2b2+s3b)].

(2.24)
Since, for our purpose, we need to compute W(+,-) and
W(+,-,-), we consider the case when there are down-crossings
in W. Suppose the kth zero-crossing is down-crossing. Then
in calculating the corresponding probability density W we
need only to take the range of integration of gi from -m to
0 instead of 0 to m. Equivalently, we may simply reverse the
sign of the (n+k)th row and column of M-1, and hence the kth
row and column of (hij) and (vi ). Thus by changing sign of
v12. in Equation (2.23) we obtain (Longuet-Higgins)

W(+,-)

=(2)-2(R11R22-R12R21)-(h1h22 -v12-12cos- v12],

(2.25)
and also by changing corresponding signs of v.. in s and b
in Equation (2.23), we get

W(+,-,-)=(2n)-3D-l (hh22h33) [I(v) +sb +
-, 11 22 33 j~ U 1b


(2.26)








Zero-crossing Interval Density Function for a Bandlimited
White Gaussian Noise


In this study the noise

white Gaussian with'a power


is assumed to be bandlimited

spectral density


for |wl 1,

for jwl>1,


(2.27)


and autocorrelation function


R(f)=f-lsinr.


(2.28)


The power spectral density and the autocorrelation function

of this noise have the shapes as shown in Figure 1.


S(w)


-1 1


R(r)

1.0




-21 -ii r


Figure 1.



The zero-crossing interval density function of the

noise is computed according to Equation (1.8)


p(0;r)=W(+,-)/W(+)-


S dt2[W(+,-,-)/ (+)]
t

and plotted in Figure 2. Experimental data by Blotekjaer(6)


S(w)= T T,

0,








21











I-








.......... ......... ... .. ...... ... .. . .............. . .. .. . ......... ........ ......




S -0






: 0 : : : :
.. .. ..... .. .. .... F ...... .. .... ... . .... ....... .. ....... ...-..








.......... ........... ... ....... ..... ...... .......... ...... ......... r......... r
E a- -
o x
o




















. . .... .... ...... .
.......... .................... ...... ...... .. ... ...... .. ......... . .. ....... .









......... .......... .......... ......... ..........











00 0 0
,
- . - - -







are also included for a comparison. The programs used in the

computation are listed under LONGO, LONG1, and MATIN. For
given 7, W(+) and W(+,-) are computed by LONGO according to
the Equations (2.20) and (2.25) respectively with the
autocorrelation function of noise

R(T)= -lsinr.

The program LONGO calls the subroutine LONG1 to compute the
term

S dt2[W(+,-,-)/W(+)],
tl
where r=t3-t.l
The computation repeats with different value of-r to cover
the range 0 LONG1 to compute the determinants of matrices which are
needed in evaluating Equations (2.17), (2.18) and I(vij)l
in Equation (2.26).








Evaluation of W for Sine Wave Plus Noise


For the zero-crossing interval density function of sine

wave plus noise, we will use the approximation


p(OT)=W(+, -)/W(+). (3.1)

The assumption used here is that the signal-to-noise ratio

is fairly high so that the zero-crossing intervals are mainly

determined by the signal frequency, and the probability of

having two down-crossings in the intervals is negligible.
Let us denote the signal plus noise by f(t). Thus

f(t)=s(t)+n(t), (3.2)

where n(t) is stationary Gaussian random noise with an

autocorrelation function of

R(r)=7-lsinr,

and where
s(t)=acos(qt+G0), with

a=ratio of sine wave amplitude to rms noise,
q=radian frequency of the signal, and

6 =initial phase angle which is unknown.

The signal s(t) is assumed to be statistically independent
to the noise n(t).
To compute W(+) we consider, as before, the joint
probability that f(t) should pass zero with a slope gl in

the time interval (t,t 1+dtI), that is, the probability of







f(t)=s(t)+n(t)=O,
f'(t)=st)'nt)+'(t)=gl (tl f%(t)=s'(t)+n'(t)=gl' 1 1

or, equivalently,

n(t)=-s(t),
n'(t)=gi-s-(t), (tl
Since the signal and noise are statistically independent,
the probability that the conditions in Equation (3.3) are
satisfied is determined by the statistics of the noise alone;
that is, it is equal to the probability of the noise satis-
fying the conditions in Equation (3.4). Thus

p(fl=0,fl'=gl)=p(n,=-s1,nl'=gl-si')

=(2y)-1|Mr 'exp(- XTM-1X), (3.5)
where
X= -s ] = -acose

g1-s1 1g+qasine ,

M= fR11 R11' i 1 0

R11' ,-R 11 ,-RO '

9=the instantaneous phase angle of the sine wave given
by e=qt+-e,
R=the autocorrelation function of the noise.

Since 60 is unknown, we may with a sufficiently large
sampling time consider 9 as a random variable which can take
all values from -v to v with equally probability. Hence from
Equation (2.2)







W(+)= (2Tr)-1 feldgiglp(nl,nl,)


=(2T)-1 r deSldglgl(2T ) j expC-iCa2cos2 +
-IT 0
(-RO')- (g1+qasine)2]]


-=(2)-2(-R0")- Jdeexp(-4a2cos2e)0dglg0
-nr 0
exp[-i(-RO")-1 (g1+qasine)23.

By a change of variable

g=(-R0") (g1+qasine),

w(+)=(2 1 )-2(_RO, )- 1 2 11
W(+)=(2T)-2(-Ro J'deexp(-a2cos2 )fdg(-R0" ) (-R0") g-

qasine]exp(-ig2),
where

c=(-RO") )qasine.

1
Letting b=(-RO") -qa,

W(+)=(2T)-3/2(-R ")I Jdeexp(-a2cos2)[
-IT
(2T)- *f dggexp(-2g2)-bsin9(2T) f dgexp(--g2)
bsine bsin9

=(2T)-3/2(-RO") deexp(-ia2 os2 )[(2T )- exp(--b2sin26)
-TI
1 bsine
-bsinO[-(2T)-_" S dgexp(-.g2)].
0
i x
Since erf(x)=2n- .fdtexp(-t2) and the average of sine from --w
0
to Tr is zero,


-Td





2 2 2 2
(2t)-exp(-ib2sin2 )+-2bbsinerf(2 in)]. (3.6)

This result was given by Rice(38).
To compute W(+,-) we need to consider the instance when

f(t)=s(t)+n(t)=O,
f'(t)=s'(t)+n'(t)=gi, (ti
or
n(t)=-s(t),
gs(t n'(t)=gi-s'(t),

The probability density for this instance is


p(nl,n2,nI'n2)=(2w)-2 -exp(-X M1), (3.9)

where
X= -s = -acose
-s2 -acos(e+qr)

gl-S1 gl+qasine

g2-s21 g2+qasin(O+qT)


M= R11, R12, Rl1, R{2 = 1, R 0 R'
R21, R22, R~1, R2 R 1 R', 0
-R11.-R .-R2 -R 2 0 ,-R',-RI,-R"

-R2I R -R1" -R 2 -R', -R";-R "
L 21' 221 -21' 22 0

Again we consider 0 as a random variable distributed
evenly from -w to w. And recognizing that the range of g2 is
from -- to 0, we write from Equation (2.2)
71 1' 0
W(+,-)=(2n)-1 rdej'dg1j dg2g (-g2)p(n" ,n2,n ,n2')
-Tf 0 -CI





27

-(2n)-3MI- JdeOdgljdg2glg exp(-iTU M-'), (3.10)
-T .0 0
where
U= -acose
-acos(9+qr)
gl+qasinI

L-g2+qasin(8+qq )

The exponent in Equation (3.10) may be written as

M' -1
-iU'M U
=-(21MI) -12a2M11[cos2e+cos2(e+q-r)]+2a2M14cosecos(e+qT)

-2aM12[ (gl+qasine )cose+[g2-qasin(e+qr)]cos(e+qr)]

+2aM13 [g-qasin(o+q-)]cose+(gl-qasin )cos(++qsr)]

+M22[ (91+qasine)2+[g2-qasin.(e +q-)]2]

-2M2 (g1+qasine) g2-qasin(9+q-r)] (3.11)

where

M11=[(-R0")2-(-R")2]-(-R0")(R')2,

M12=R'-(-R" )R+(-R")],

M13=-R'[(-RO")-R(-R")]+(R')3,

M14=-R[(-RO,)2_(-R")2]+(-R")(R')2

M22= (-R")(1-R2)-(R')2

M23=-(-R")(1-R2)+R(R')2,

!MI =[ (1-R)-RRo")+(-R" )]-(R' )2][ (I+R)[ (-RO")-(-R")]-(R')2





28

We make a change of variables

(2 MI )-gl=rcos(X+~),

(2)MI) -g2=rsin(kx+I),

wnich'yields, after some rearrangement of terms,

w(+,-).
= 0/5 -3M3 V SIT 3 2
-=-31M3/2 dofdr 'dkxr cos(2%)exp[-(Dir2+2D2r+D )],
-T 0 --T
(3.12)
where
D =M22-M23cos(2x)

=(M22-M23)cos2 +(M22+M23 )sin2 X

D2=alMI ([-[ (M12-M13 )cos(qr)+q(M22-M23) sin (-iq)])

cos(e+F-qr)cosh

+[ (M12+M 13 )sin(iqr) -q(M22 V 23)c os (qT) ]sin( 8+2qT) sink],

D^=a22MI -1[[(M11+ 1,4)cos2 (q-)+2q(M12-M13)sin(iqT)cos(GqT)

+2 (M22-M23) sin2 -r) ]cos2 (+1qr)

+[(M11-M14)sin2 (-q-) -2q (M12+M13)sin(qr)cos (qT)

+q2(M22+M23)cos2( q-) ]sin2 (e+2q)].


We simplify the above expressions by writing

C0=R',
C,=1-R,

C2=1-+R,
C =-RO"+R",





29


C4=-R0 "-R",
C=(I-R)(-R-R")-(R')2=01C 4-002,
C6=(1+R)(-R"+R")-(R' )2 =C2C-C2,
=C7COcs(iqT)+qC2sin(iqr),
C8=C0sin(iqT)+qC1'cos(i-qT),

and replacing (e+qr) by e to obtain

D1=C 25 os2X+C1C6sin2 ,

D2=-a(C5C6)-'(C5C7cosecosx+C6C8sinesink),

D3=a (CC6 )-1[C2 C5(C 72+Cgcos2( 2qT)
+C 1 2 2(lq-) )in 0]
+C1-1 (C82+C5sin2(iqi))sin2,

IMI=c5c6.

By integrating the Equation (3.12) with respect to r
we finally write (Cobb)
W(+,_)=(2y)-31MI3/2 '2 r 2
W(+,-)=(2)3 M/2 de dxD1-2cos(2X)[(l+D 2)exp(-D3)
-I -r 2
-nD4 (3/2+D42)(1-erfD )exp(-D5 ),

(3.13)
where
D4=D2/D1 ,


D=D 3-D 2







Zero-crossing Interval Density Function for
Sine Wave Plus Noise


The zero-crossing interval density function for sine

wave plus noise is computed according to the Equation (3.1)

p(O -T)=W(+,-)/W(+)
with the same noise characteristics as before. Several of

these density functions are plotted in Figures 3-8. STOSN,

SINO, and COREK are the programs used to compute p(0;r).

The computation of W(+,-) is carried out by SINO according
to the Equation (3.13). The program COREK computes Equation

(3.6) to evaluate W(+) for the signal-plus-noise case. The
data computed by SINO were stored in the computer operating

system by the program STOSN, and later retrieved by COREK

to complete the computation.
The second term in Equation (1.8)

S dt2[W(+,-,-)/W(+)] (3.14)
tl
represents the probability density of there being at least

one down-crossing between an up-crossing and a down-crossing
separated by an interval T=t -t,. When the amplitude of the

signal is high, the zero-crossing intervals of the signal
plus noise are close to the half period of the signal wave,
Th=n/q, and the probability that the zero-crossing intervals
are much longer than Th is very low. In other words, if r is

larger than Th, the probability of having at least one down-

crossing within the interval r is very high.
The zero-crossing interval density function in Equation





31.


(3.1) therefore contains, with the omission of the second
term in Equation (1.8), erroneous distribution nearr=3Th.

When the signal radian frequency q is 0.6 or higher, the

erroneous distribution falls within the interval (0,15) of r.

These erroneous data are discarded, and the main distribution
curves are extended smoothly to vanish at r=15 by the program

DPOL.

The zero-crossing interval density function (3.1) is a
good approximation if the ratio of sine wave amplitude to

rms noise is large. In the cases when the ratio a=l or 2, the

error in the approximation is fairly high. Instead of compu-

ting the density functions for these cases directly from the

Equation (3.1), we interpolated the data by the program

SAPOL from the density functions for a=0,3,4 and 5 using a

spline interpolating function (Appendix D). The zero-crossing

interval density function for a=0 is just the density function

of noise alone.

The program INPOL expands the data points of the density
function by locally interpolating with a spline function. The

area under the computed zero-crossing interval density func-
tion usually has a slight deviation (5%) from 1, and this

is corrected by normalization.






34






*
.......... ....... ... ..... ....... ...


....... ................ ........ ........ ........ ........ ........ ....... ........























0 -; C a
II n : : "





...... .. .. ... .. ... .., .:















a. .

,-




r1\
0)












a
". ....... .... . . .... :........ ....... ....... ........ ........ . ..... .



........ ....... ........ ... ... ....... ..... .... .. . -



S'r
....... ....... ........ ........ ........ ....... ........ ....... !. .........V.. o
v
















0.
S......... .. . .

o o o d o o o o 0
c. d d. .--







35

.. ....... ...... ....... ...... ... .... ........ ........ ....



at







....... ..... .. ....... ...... ......... . . .




U-g
o o-




















So B ao t- a V C
o



























0


t- P











II 11


i... + . ....... ..... .. .... ... ....... ... ..... ... ...... .....





















a o o j n .- V r-

0.


~II_~






S36


....... ............... ....... .......................................... ...




...i
a-

'..



o cr
....... ........ ........ ; ....... .. ........ .... ... ........ .... ... ....... . .























,o ,o o ,d d o o ,
a. ..





-I
....... ..... ... .. .................... ............... ....... ....... ......... ........ .



























*...... ......... . . .. ... .
CL











..V 0 9*0 ... ... .........


S0 0 0










C
..-.... ..... ........ ....... .. ..... 0...... 0 0
a-
. . . ; i ; . . . . . . . . . . . . ^ . . . . . T


V ^














O ~ f' o 0" S 0 0 0 0
a.-







37


__ ..... ....... ... ........ ,..................... ........ ...........



r-
.. ..... ....... . .









* *. .. .. *...






............. ... ..... ... .
oto













.......... ........ ........ ..... ........ ................ ........ .. ........



















r.
'" O i











..a. . . . . .


























........................ .... ...... .. ......... 1
0 0


























, c to o o o o o o o
o.0
co* U oD














CHAPTER III
IMPLEMENTATION AND EVALUATION OF THE DETECTOR



Implementation of the Detector



In implementing the detector it is assumed that the

signal-to-noise ratio and the length of a signal are known

to us. We also assume that the signal-to-noise ratio is fairly

large so that the additive noise on a signal only disturbs

the zero-crossing interval determined by the signal frequency

and does not add more zero-crossings in the process. Then we

know the number of zero-crossings needed to observe the signal

in full length.

The schematic of the detector is shown in Figure 9.

Zero-crossipgs of f(t), both up-crossings and down-crossings,

are first detected and the output pulse of the zero-crossing

detector is used to strobe the system. The intervals between

zero-crossings are measured and compared to predetermined

thresholds. The output Zi of this comparator has values such

as

Zi= 1, if T
S0, otherwise, (4.1)

where

T0=the lower limit of the threshold of the zero-crossing

interval comparator,






39











0



PI 0
po

c

O

*-I 0 0













I C

0 *
C)l 1 0

011




*4







m00C
I Cd


















k Oz>
o "-I








T1=the higher limit of the threshold of the zero-crossing

interval comparator.

The outputs of the zero-crossing interval comparator are

stored in a shift register and the content of the shift regis-

ter is counted by an up-down counter. If the signal s(t) is

composed of n cycles of sinusoidal waves, then there are 2n

zero-crossings generated by the signal. To observe the signal

in full length, the counter sums up the content of the shift

register in such a way that

0, if i=O,

C.= ci-' if i-2n' (4.2)
ci= (4.2)
1 Ci -+1, if Z =1 and Zi-2n=0

Ci_1-l, if Zi=O and Z =1,

where Z.=0 for i(O. Therefore the value of Ci at any moment is

i
Ci.= -2 Zk (4.3)
1 k=i-2n


and this is compared to a threshold to determine the presence

of a signal. Thus if we set the threshold to m, then the

output of the comparator is


Di= 1, if m
0, otherwise. (4.4)

Let us denote that

H1=the hypothesis that a signal is present,

HO=the hypothesis that a signal is absent,

and for given a and q of the signal,







Ps=P(Zi=11H1)=the probability that when a signal is
present T. lies within the interval (TOT),

Q =P(Zi=0|H1)=the probability that when a signal is
present T. lies outside the interval (TO,Ti),

Pn=P(Zi=1lH )=the probability that when a signal is
absent Ti lies within the interval (TO,T1),

Qn=P(Zi=0IH0)=the probability that when a signal is
absent Ti lies outside the interval (TO,T1).
Then
Ps -s=l,
P n+Q 1,
and Ps corresponds to the area under the curve bounded by
T=T0 and T1 in Figure 10(b), and Pn corresponds to the area
under the curve bounded by the same limits in Figure 10(a).
To simplify the analysis, the detection probability of
a signal is approximated by the detection probability deter-
mined by comparing threshold at the end of the signal. This is
the lower bound of the detection probability and is given by

Pd=P(Di=1 H )=P(Ci>m| H)


2n 2n k 2n-k (5)
k=m s s

The corresponding false alarm rate is

Pf=P(D =1 H )=P(C im|H )

2n 2n k2n
= m ( )Pk Q 2 (4.6)
km n n






42

p(o; r)



0 .4 ........... ........ ........... ......... .. ................... ............ ... .... .






0"1" noise olo ie

0. 6 10 1 14






T0 T,

(a),











0.3 ... ... ....... .... ...... ... .......... .......... I ..........,
0.1... .


















2 4 6 B 10 12 14 16
To T, T

< b)






Figure 10.
Fi~ur-e 10.








Now the objective is to maximize the detection proba-

bility Pd for a given false alarm probability Pf. This
strategy is known as the Neyman-Pearson criterion of decision
making, which does not involve either a priori probability
or cost estimates (Whalen).
Let r=Ps/P Then
S2n

d ( )(rP )k(-rP )2n-k (4.7)
k=m

The binomial expansion


f(n,m;x)= E (n)xk(-x)-k (4.8)
k=m

is related to the incomplete beta function by


f(n,m;x)=[ Sdtt-(l-t)n-m]/[ dttm-1(l-t)n-m (4.9)
0 0

(Abramowitz). Since the right side of the above equation

increases monotonically with respect to x for given values
of n and m, the binomial expansion on the left side also
monctonically increases with respect to x. Then we see that,
for a signal with parameters a, q and n, Pd in Equation.(4.7)
is maximized by maximizing the ratio r for given m and Pn'
and thus Pf. The maximization of r is done by the program

PCOM3.
PCOM3 integrates the zero-crossing interval density
function of noise from T=0 until the integrated area equals

Pn, and the zero-crossing interval density function of signal








plus noise is also integrated for the same interval to obtain

the corresponding P This integration procedure is repeated

with an incremented lower limit of integration until the
range of t is exhausted. Among the Ps's obtained by such

integration, the largest value is picked out as the maxi-

mized Ps for a given P The integration interval correspond-

ing to the maximized Ps is then the optimum threshold for

zero-crossing interval as in Equation (4.1).

Figures 11 through 14 show the relations between P and
the maximized Ps for different values of parameters. In the

Figures 13 and 14, P is plotted against q with fixed para-

meters Ps and a. From these graphs we observe that q=0.4 in

general gives a large value to the ratio r=Ps/Pn.

The final step of maximizing the detection probability

Pd for given Pf is the determination of the threshold m. The
program RCOM computes Pd and Pf for different values of P

and m according to the Equations (4.7) and (4.6) respectively

as shown in Figures 15-19. From these graphs the optimum

value for m could be determined to give the maximum Pd for

given Pf of a signal with given parameters. For instance,
the optimum value of m for a signal with parameters a=2,

q=0.4 and n=4 to be detected with maximum probability for

a given false alarm probability Pf=10-2 is 6, while 7 is the
optimum for Pf=10 .







45

U ..
..... ............................... ...
















.O






















0o 0 r







a-1
S. ...... ..

oo



















or











0 j . .. .. . ....
o 0 0 0 0 0 0 0 0
o o o
























0 0 /-



.............. ............................... .......... .. .. 0








0j




















a a

: a









.................. ........ ..... .




. ..... ........ ................... ........... .
II II












.. . *
... m






o -- a e r- a v 4 ,
w' .i. .. a; ; o o n-1 ~






33

. ... ...... ..... ........ .......... ... ........ ....... .. ... .... ...












........ ........ ...... ...... . ... ..... ..... .



............. . ... .. ........ ....... ........ ........ ........ ...... ...... .







.......... ........ ........ ........ ........ ....... ........ ...... ....... ....... ..










:o o o o o o
'.-4















*r4
a-e


































....*. . . . .. ......... . ............ ....... .... ... .



............ ........ ................. ........ .......... ........ ....... ....... 0



................ . ........ 1 ........ ........ ........ ...... ........ ... .... ..
'Cr --I-rT--I-i--i-- i--1--i--r-*
- o m I I I?
o CCoo N cooo -
4.



























I :-
'.4



















-I . ~.......,.......i....... ..




1- - r rr-r i r-r

0 S fl Y f .
.5 a a C a 0 0
a-cr o n








46

c
a-
......,..... ....... ....... ................. .






... ............... ........ .. ..... ........ ....... ....... ........
o



........ ............... ........ ....... ....... ........ ....... .
0








.. ..... ..... .. ............. . .................


................. ................ ....... ......... ...... .. ........ .....
............... ......... ......... ................ .......... ................. ........





o:0


















C(




-
0

I Cf
"""""""~""""""""" ""~" ~a-


- 0 0 0 0 0 0 0 0 0
Ocrec-Drwnnr.
000000000o


































I I I I I I I I I I

* 0 0 0 0 0 0o p



Co

f d





.... .. .
........ ....... .................... : .... ....... ................ ................









'n :
o.i




o o o
o. o..o........ .. ................. o










\ 0 0 0 0 d 0
J ..... I .... ... .... .L .. ,....


.. .... ./ :... / ...../ ... .. ., -.:

i i 1 r T i i i i -








S48



O"
S....... .... ... ....... ............. ... ............. ... .... ....... ....... .








. .* . :











............. ... ........ .......... : ....... .... .. .. .. .... .. ....







..... .. ........ ................. ........ ........ ..........




I* C-
0 0) 1 0 0
C.













*d




.. .. .... '.....................;.... ...............








....... ........ ....... ......... .... .. ... .. .. .... .. ... .
..... . .. ... . a








































it o






.. ... .... o.............
a. a0

C)k

0r


0 C I r C In V ) N

0 0 0 0 0 0 0 t-.
a;- ..... .. -..



m0
a~V~.s o j*0
j...~..i. i -;in..


















0O
))- .....!......)......:...C)


Si -
-r- ----r-r-- --r--
-* 0 0 0 0 0 0 0 0 0































.. -E.
........ ............... ........ ....... ...... ................. ........ .......









_4 .............


.............. ... . ........ . . .. .
0 ,


II iI n








































...... :. ............... ...- ....... ........ ................ ................. o
o o0



































.. .. .... ... ...... ..... ..............jE .... . .. .............. o
: : : 9' ^


r --
0 0 0 ~i- 0N )
0 0 0 0 C) a a








50










.. .. ... ............... ........ ........ ....... ..... ............... ....;. o









........ ........ ....... ... .



...... .... .. .... 0





. . ...... ...... ..
I - i- -i -- i -- i -- i -- ; -- i i --















.o .. . ......... .
0
I





*1 * N ** *
0 0 0 0 0 0*
oitrcl i '0
o-

L


r"



O~n(D CY I~O0

d o~ d d d d0

Uc


I I I 1 I I I I I I
o0 Ct C C o o o N s -
- 0 o 0 0 0 0 0 0








51


.. ...... .............. .......... ............. ..... ... .......




.-.... ...... ...... .... ........... ... ....... o


























o o
0 f I> .Ii
-0















T-' o o ,.o o o o 6 o o
Ii i a

... 0






oImII l s








*r



........ ......... ...... ............. ........... -
-o

























........................ .... ..... ............ ....... .. .
.. -... ..... ........ ....... 0
E .


......~... ......... .. o






.... ... ..... ...... ........... ... ... ..... ... .. ..... .... ...... .. ....... .
a*


























........ ....... ........ ......................... ........ ........ ........ ........ 10




.o. ... ... .... . . .. .
00* 0 0 :
0


0m 0 0
fl-a







52


?o



.. .... .. ... ....... 0






. ....... .... .. ..... .......... ....... .- ........ ........ 'o
.. .... ..... .. *. .......... ... .--...*











"" 1


V: 0 co0

o o" c .-
. R... . ............." ....... ........ .. . . . . . . o







at CO F! c!
r*r













1 c





0








.. ....... .... ...... ........ . ..... .. ........ ........ ... . ... -. 'o
o. "
:


i- I i ... I 'T--"-
o a e i *
- 0 C o0 0 0 C 0








53

I f
O


o














................ .... ... .... ........... .... ...... ... ............. o







.......... ... .... ..... ..o



H1 II I o
: c:







,0 :0 C :





















....... ................ ........ ......... ................ ......... ................. o


................ .. .. . .

..... .......... ..............................
S. .. o.


010








II :if oli
2* r
o




.. .. .... .... ... ............ ........ ,.. .. ...... .... .. o





-a>
a

II ii II

00-C *


rr~ rrnmr

0 0 C V C. c
0 0 ~ 0 0








Comparison to Correlation Detector


The performance of the detector is compared to the

conventional correlation detector. Suppose that the para-

meters of the signal s(t) are completely known to us, and

let H1 and H0 denote the hypotheses that the signal is

present and absent respectively. The correlation detector

chooses Hi if

T
Jf(t)s(t)dtV, (5.1)
0

where the interval (0,T) is the duration in which the signals

are received and V is a threshold. The implementation of the

correlation detector is shown in Figure 20.





T H
f(t) X --- Threshold or
0 H1

s(t)


Figure 20.




A correlation detector is equivalent to a matched filter

followed by a threshold comparator (Whalen). This equivalent

detector which implements the operation of the correlation








detector using a matched filter with impulse response h(t)

is shown in Figure 21.






f(t) h(t)=s(T-t) Threshold H or
0 t t=T




Figure 21.



The decision rule in Equation (5.1) stands on the

assumption that the noise is white Gaussian. The threshold

V is determined to meet an employed detection criterion for
the particular application. To make the environment similar
to our previous analysis in analyzing the correlation detec-

tor, we assume that the noise is bandlimited white Gaussian

with power spectral density

S(w)= ) for Iwicl,

0, otherwise (5.2)

and autocorrelation function

-1
R(T)=-r1sinr. (5.3)

The Neyman-Pearson criterion, which involves neither a priori

probabilities nor cost estimates, is also employed here. The
strategy of this criterion is to maximize the detection






56

probability for a given false alarm probability, and this
could be accomplished by using a likelihood ratio test
(Whalen). We approach the analysis with the case where the
signal is sampled at m discrete times.
Let us denote the m sampled values by

fk=sik+nk (k=l,...,m) (i=0,1), (5.4)
where

0k=0, Slk=s(tk), nk=n(tk), fk=f(tk)
and tk (k=l,...,m) are sampling times.
If a noise is sampled at m discrete points, the
covariance matrix M=(m i) of the noise has the elements


mij=E[ni.n]=R(t.-ti) (i,j=l,...,m). (5.5)


Since the noise is Gaussian, the joint probability density
of the m sampled values of noise is


p(n)=p(n,n2,,...,n )

=(2u)-mlMI exp(-X TM-1X), (5.6)

where X is the column vector of sampled values. By writing

Q=M-1.(qi .)- (i,j=1, ... ,m),

the Equation (5.6) becomes
1 m m
p(n)e=(2)-I MI exp(- S qnijn n). (5.7)
i,l=1

Since nk=f.k-shk (h=0,l), the above equation could be







written as

p(n)=p(f-sh)

=(2T)-mlMJj expL[- q(f -sihi jf-shj
i,j=1 h h j

Then, when the signal is absent in f(t), h=0 and the joint
probability density of m sampled values of (fk-sk)
(k=l,...,m) is

PO(f)=p(f-s0)=p(n)
1 1
=(2n)-MIM-exp[- E qij(fi-oi)(fj-oj)].
i,j=1
(5.8)
If the signal is present in f(t), then h=l and the
corresponding joint probability density of the m samples

(fkslk) (k=l,...,m) is

Pl(f)=p(f-sl)=p(n)
1 1 2
=(2)-_ImlJMI -exp[-i qi i-sli )j-f- s
(5.9)
The likelihood function for the hypothesis H0 is
the joint probability density function p0(f)=p0(f1.....,f )
Likewise pl(fl,...,f ) represents the likelihood function
for the hypothesis H1. In terms of these likelihood functions,
a likelihood ratio r is defined by

r=p1(fl ..., fm)/po(fl...fm)
m
=exp[-i,= i q.(f -sli)(f-slj)]/

m
exp[-i E qi.(f -soi)(f -s
i, =1 J a







m
=exp[-j E qj(2f s0j-2f slj-sOi1s+s5lisi), (5.10)
i,j=1

and the likelihood ratio test is to.choose

H1, if r>,ro
HO, if r
where r0 is the decision threshold.
If the noise is not white, the evaluation of the likeli-
hood ratio r in Equation (5.10) involves a matrix inversion
from the covariance matrix M. When the signal is sampled at
a large number of points, this matrix inversion is not
practical in considerations of speed and simplicity of the
system. This is also the case with our bandlimited white
noise, and we simplify the analysis by sampling f(t) at the
time intervals for which the sampled noises are uncorrelated.
Since the autocorrelation function of the noise has
zeros at T=kT (k=-1,2,...), if we sample f(t) at the inter-
vals 6t=T, the samples are uncorrelated; that is

E[n.n ]=R(t j-t)= j R(O), if i=j
0, otherwise. (5.12)

Then, since, from Equation (5.3), R(0)=1, the elements of
the matrix M in Equation (5.5) become

mi =E[nn ]= 1, if i=j,
0, otherwise, (5.13)
and thus,
M=I,






59

and

qij= 1, if i=j,
0, otherwise. (5.14)

Then the likelihood functions in Equation (5.8) and

(5.9) simply become
i m 2
po(f)=(2)-imZexp[- (k0k)2 (5.15)
k=1

and


pl(f)=(2)- mexp[- (k-slk)2]. (5.16)
.k=l

The likelihood ratio r in this case is

r(f)=pl(f)/po(f)
m m
=exp[-i Z (f k-sk)2/exp[-i 2 (f-Sok) J (5.17)
k=1 k=-

or, upon combining terms,

r(f)=exp[-2 E (2fks0k-2fks k-s0k +sk2)]. (5.18)
k=1

The decision rule is to choose H1 if rr0, or equivalently,

by taking natural logarithm, choose H1 if
m 2 2
E C(kslk-f ksOk+sk -~slk )n(ro). (5.19)
k=l1

Let
m
G= E (f ks-kfksOk Ok2-s lk2). (5.20)
k=J.
To evaluate the performance of this detector we need to know
the distribution of G under each hypothesis. When H0 is true,






60

fk=nk and the expectation of G is
m 1 2
EO[G]=E[ E (nkS1 -nkS0k+S0kj-i lk )]
k=l

Since Sok=O,
m
Eo[G]=-1 8sk2
k=1

The variance of G under H0 is

Vo[G]=E[(G-Eo[G])2].

But

G-EO[G]= E nkSlk'
k=l

so
m m
VoEG]=E[ Z Z nknj s!kslj
k=l j=1

m m
= E E[nkjn.]slkslj.
k=l j=l,

Since

E[nkn j=' R(O)=I, if k=j,

0, otherwise,


it becomes
m
VO[G]= E s lk
k=l

Define
m
S= s lk2.
k=1

Then,

Eo[G]=-iS,


(5.21)


(5.22)


(5.23)


(5.24)



(5.25)








VOl[G=S.

Therefore the density function of G is


p0(G)=(2Tw) -S-exp[-(G+iS)2/S]. (5.26)

In a similar fashion we obtain the density function of

G under the hypothesis H1. If the hypothesis H1 is true,

then fk=slk+nk and

m m
E1[G]=EL[ (Slk+nk)slk E (-sk2
k=1 k=1

m 2 m 2
= lk -i k s1k
k=l k=1

m 2
= s (5.27)
k=1

m m m
G-E1[G]= E (slk+nk)slk+ (-s )- 12
k=1 k=l k=1
m
= E nksik.
k=1


VI[G]=E[(G-E1[G])2]

m m
=E[ E Z nkn slkSlj]
k=l j=1

m.., .m
= Z j E[nknj]slk sl
k=l j=l

m 2
Sk s (5.28)
k=1

Thus


p,(G)(2- 2S 2exp(L-(G-iS)2/SI].


(5.29)








Since we choose H1 if G.ln(ro), by writing c=ln(r0),

the false alarm probability is
00
Pf =P(D1 HO)= FdGp0(G)


=(2T) -S~_- JdGexp[-i(G+iS)2/S].
0


By a change of variable z=S~2(G+IS),

Pf =(2Tr) f S dzS2exp(-z2)
(c+-S)/S2


=(2)"- J' dzexp(-z2). (5.30)
(c+S)/S

The probability of detection is
CO
Pd=P(D1 H1)= 1 dGpi(G)
-C

00 2/s i
=(2)"-S-' JSdGexp[-2 (G-S) 2/]
c

1 0 1 2
=(2n)-S- ~J' dzS exp(-iz2)
(c-IS)/S2


=(24)- dzexp(-iz2) (5.31)
1
(c-IS)/s'

The performance of the correlation detector is computed

by the program PLIN and shown in Figure 22. We note that the

performance of the correlation detector is superior to the

zero-crossing measuring detector in the assumed circumstances





63

of the signal and background noise. For instance when a=2,

q=0.4 and n=4, the false alarm probability of the correla-

tion detector at Pd=0.9 is 4x10-7 while the false alarm

probability of the zero-crossing measuring detector at the

same detection rate is 3x10-3.

This superior performance of the correlation detector

to the zero-crossing measuring detector is based on the

following assumptions about the signal and noise: the

frequency and amplitude of the signal are completely known

and fixed; the noise is stationary in Gaussian distribution.

If the practical situation deviates from the assumptions,

then the performances of the two detectors will also deviate

from the previous results. The merits and operating charac-

teristics of the two detectors are further discussed in the

next chapter under such practical considerations in detecting

signals from sleep electroencephalograms(EEGs).










-ocC
.




..... .. ........ ......... ........ ...... ... ..... ........ 0
. ............ ................ ........* * .O




















0 0


............;.. . . -





. f O f. .. . . ,,. .1 .
S II o o o .

0 0 0 0 0 0 0 0










*r






.. -.. ..: ......... -




... ... . ..... ....... ..... .. -........ ....... ......... o






















j-" -- ----,----T----T-- --'--,F---j---

o a 0 t ,- o .m V C V -,
. 0 o0 0 0 0 0
0:













CHAPTER IV
FURTHER DISCUSSIONS AND CONCLUSION


Operating Characteristics in Practical Situation


The zero-crossing measuring detector has been applied

to the detection of phasic events in sleep electroencephalo-

grams(EEGs). Examples of EEG activities from a cat are shown

in Figure 23. The arrows indicate phasic events in a sleep

EEG. Among various phasic events, alpha, sigma, and beta

waves could be modeled as signals of sine wave form imbedded

in the noise of background EEG.

The physiological origins of the phasic events as well

as the irregular activities of EEG have not been completely

understood as yet and neither the statistical characteristics

between the phasic events nor the irregular activities have

been studied. Even if we can assume for practical purposes

that the phasic events are independent processes superimposed

on the background EEG, the nonstationarity and the lack of

Gaussian property in the sleep EEG activities make it diffi-

cult to apply the previous analysis to the quantitative

determination of the advantages of the two methods in detect-

ing phasic events from sleep EEG. Hence in the following we

attempt only to make qualitative observations on the merits

of the zero-crossing measuring detector.






S66


*b -* 1 *








-.- --.-
_______________ -S -
4-5




_---- _- _ _--- _--- -- ^

- 4- --?:yS-s-r
r

~~-v

-- .. ----

-ii
--- -- r-- -
-'-^ _- ,









.v 4 ~ *, ^ __.._ _. _= .,__ ...__
*'_ls='ji3 -^?-





-_ --- -- = ""
-4--. --- s;- --
--^------ --~-. --^ -....-~
-^S" _L ~ ___ "''Sp _; ___ __ ^^ __ "^ --- i __ -""








_?~~ *~ '' e"".- f.





<1^ --- ^ *-* S 4^ *- **
-> ^. r ^ . ^ -








When we attempt to detect phasic events from a sleep

EEG we encounter the the following situations: 1). The

frequency of the phasic events varies among individuals. The

frequency of alpha activity, for instance, ranges between 8

and 12 Hz, sigma activity between 12 and 15, and the a priori

probability density of the frequency distribution is unknown.

2). The amplitude of phasic events varies due to the elec-

trode resistance and individual differences in age and

physical parameters. 3). Background EEG activity level varies

and artifacts are introduced by movements. The background

EEG activity varies also due to the differences between

individuals and, within an individual, the variation of

activity inherent to a sleep process. Examples of movement

artifacts recorded from a cat are shown in Figure 26.

To qualitatively observe the operating characteristics

of the two detectors in the described situations, we do the

following analyses:

1). To see the effect of the varied frequency of the

signal on the performance, we compute the detection proba-

bilities of the detectors to signals of different frequency.

Suppose the correlation detector is tuned to signal s1, and

we receive a signal s2 with a same amplitude but different

frequency from s1. Then, since fk=s2k+nk, G in Equation

(5.20) becomes
m m
G= E (sk+nk)sk E s (6.1)
k=1 k=l

The expectation of G is







m m
E2EGJ= E s2kslk- 2 slk2. (6.2)
k=1 k=1

The variance of G is


mm m
V2[G]=E[(G-E2[G])2]=E[ E E nkn slkslj
k=1 j=1

m slk2=S. (6.3)
k=1
Thus the density function for G is


P2(G)=(2Tr)-S-expC-(G-E2) /S], (6.4)

where E2=E2EG].

Since the detector is tuned to the signal sl, the
threshold remains the same as for s1, and, hence, the detect-

ion probability for s2 is

00
Pd=P(Dl H2)= JdGp2(G)
c

=(2z) -"S J'dGexp[-i(G-E2 2/S]
c

=(2T)- X S dzexp(-iz2), (6.5)
1
(c-E2)/S2

where H2 is the hypothesis that signal s2 is present, and c

is the detector threshold set for signal si.

The detection probability of the correlation detector
for signals with a different frequency is computed by the

program PSEG with the following parameters:
a=2 for both sI and s2,

n=6,







c=threshold set to detect s1 with a probability of 0.95,

that is, set for P(D1jH1)=0.95,
where a is the ratio of signal amplitude to rms noise, and
n is the length of signal in number of cycles. The result is
shown in Figure 24(a). The detection probability for a signal
with a different frequency, P(DjIH2), sharply drops off as
the frequency of the signal slightly departs from the tuned
frequency of the detector.
The detection probability P(D1IH2) for the zero-crossing
measuring detector is computed in the following way. Since
the detector is tuned to sl, the thresholds m in Equation
(4.4) and TO and T1 in Equation (4.1) are already determined.
From the Figures 16(a) and (b), we note that the optimum
threshold for m is 6 for a signal with a=2 and n=4, and m=12
when a=2 and n=8. So with the present signal parameters of
a=2 and n=6, we assume that the threshold for m has been
determined at 9. As before we want P(DI H1)=0.95, which
consequently determines the optimum thresholds (T0,T1) for the
zero-crossing interval.
The probability of detecting s2 is then

Pd=P(Dj1H2)
2n
= k V 2k(P2)2n(6.6)
k=m

where P2 is the area under the zero-crossing interval density
curve of s2 plus noise bounded by the interval (TO,T ). The
computation is carried out by the program PCOM4 and RCOM,
and the result is plotted in Figure 24(b).






70



: .... ............... ... ...... .. ............ : ....... .......... ........ ................ .
0 0








... ...... ..... . ... . . ... ....

I...!
-0 . . ..... ..... . .




: to : io


.......... ....... .... ........ ....... ........ ........ . ........ -P .. .... .. ...... "




I 0 ,a 4"




acE ,.
. ......... .... .. ......... .......... ...... .. .. ..... . .... ....


* o E











0-
.... ...... i ...... ........ ....... ... .. ...... ... .... .. ... ......







S E .. ..... ........ ..... ... ...........
.......... I .... ..... .:....... .....:,... . .. .. :

".....' .... ... ..j ...... j ........ ........ ................

: o !







... ... ......... ....... ........ ........ ... ........
p : : i ; - -















T a 0
S I I I : I







We note from the two curves that the operating charac-

teristic of the zero-crossing measuring detector is more

suitable to cover the broad range of signal frequencies. To

cover the same range with the correlation detector without

affecting the detection and false alarm performance, we need

several channels of detectors in parallel, which is not very

practical.

2). The zero-crossing interval distribution is not

dependent upon the amplitude of the signal and noise so long

as the signal to noise power ratio is maintained constant,

which is the usual case with sleep EEGs. And, therefore, the

operating characteristic of the zero-crossing measuring

detector is not affected by the varied strength of signal

and noise.

However, a reduced amplitude of a signal greatly decrea-

ses the detection probability of the correlation detector.

The detection probability of the correlation detector for a

reduced signal can be computed in the same way as done for

varied frequency. In this case the first term in Equation

(6.2) becomes a function of amplitudes of two signals instead

of frequencies, and the detection probability for the reduced

signal is expressed again by the Equation (6.5) which could

be computed by the program PSEG with a few modifications.

The changes in performance of the two detectors are

shown in Figure 25(a). The horizontal coordinate is the ampli-

tude ratio of the reduced signal s2 to the tuned signal s1.

The frequency and the length of s2 are assumed to be the same







72

-c

I I*.
... ............. .. .. ........ ..... .................... .. .......... .......
I n :O
.0


i..... J I ..... ... ...... ....... .
























-oo
I 0 Ii "





o o t o C -- _
..... ... ... ........ ........ .. .. ........I ....... .















S0 0 0 0.
a-a
.....1 :......... ....l... ........ ......... ........ ..


L,0












-C0
...... ......
....... .... .. .. ... .... ........ ....... ..... .. ... ... ..... ........ j






























00 0 u0 0




*. .* '.. . ..- ..- .
.. ....... .......... ....... .'-- ........ ....... I ......... .
3 0 ,


...... ........ .... ........ .. .... ... ....... .... . .......... . .i. . .









1 I
. .. . .. . . . . .









- -"-








as s1 and the signal-to-noise ratio is assumed to be constant.

The detection probability of the correlation detector for

signals of reduced amplitude drops off rapidly as the ratio
decreases while that of the zero-crossing measuring

detector remains constant.

3). The increased noise level, assuming that the noise
characteristics remain the same, does not increase the false
alarm rate of the zero-crossing measuring detector. Also,

since the signal to noise power ratio usually remains
constant, the detection probability of the signal is not

affected.
However, if the noise level is increased in the correl-
ation detector, the false alarm rate also increases to a

considerable degree. Suppose the power level of the noise
2
has been changed by a factor of h2. Then

E[nknj]= h2, if k=j

0, otherwise, (6.7)

and, since fk=nk,
m
E[G]=- E sik2=- S, (6.8)
k=1

V[G]=E[(G-E[G])2]=E[ E Z n snlksi]
k=1 j=1 k

=h2 E s k2h2S. (6.9)
k=l

Thus


PO(G)=(2w) h1 Shexp[-I(G+kS)2/(h S)],


(6.10)








and the false alarm probability is
00
Pf=P(D1 HO)= fdGpo(G)
c


=(2T)4Ih-1S- .dGexp[-(G+AS)2/(h2)
c

=(2T)-~ dzexp(-iz2), (6.11)
(c+iS)/(hS)

where c is the threshold set for the desired detection

probability of the signal. The effect of the increased noise

level on the false alarm probability of the correlation

detector is computed by the program PNOI according to the

Equation (6.11), and plotted in Figure 25(b). The false alarm

probability of the correlation detector increases rapidly to

0.1 and gradually approaches to 0.5 as the noise level

increases. The false alarm probability of the zero-crossing

measuring detector remains constant regardless of the increas-

ed noise level.

Since the false alarm probability of the correlation

detector is susceptible to the noise level, the detection

threshold of the correlation detector could not be lowered

arbitrarily to increase the detection probability of weak

signals. Consequently, the optimum threshold of the correla-

tion detector for detecting phasic events in EEG should be

determined according to the individual strength of the EEG

activity.

On the other hand, the zero-crossing measuring detector

does not need any threshold adjustments even if the strength








of the EEG activity varies over a wide range due.to indi-

vidual differences or electrode resistance.

Typical examples of the detection performance of the

zero-crossing measuring detector are shown in Figure 27. A

square pulse is the output of the detector indicating the

presence of a signal. The immunity of the detector to high

power artifacts is demonstrated in Figure 26.



Conclusion


The performance of the zero-crossing measuring detector

in an idealized environment of signal and noise was inferior

to that of the correlation detector. But when the conditions

of the signal and noise deviated from the assumptions, the

zero-crossing measuring detector maintained the quality of

performance while the performance of the correlation detector

was susceptible to the change of environment. The sleep EEG

is a very complex process; both the signal and noise have a

wide range of variation. When the zero-crossing measuring

technique is employed to detect phasic events from sleep

EEGs, it should provide a detector which is excellent in

performance and applicable to a broad range of subjects with

a simplicity of implementation.















_-I_


*`- -






2-:'-- -E,.-





-- ----





''1--- --I. -- -



--k-. -











-- --r-

F. I~


-wc







I.






,-


















_ -I--_-










~--- -


r;eir
---~I~--~
----------"7i--e--
--Wc~-
---~a
--=~;5=-
-4~-t
-~CL ~


~:~S~L----
L ----- T;

i
~~. .~---~r~-;r_
`i: ---------~c






--L-_-~L-~












-- -------~~;~~---------

--,

--




---- -- ----- ----r~s

----- ---~-~1 ____












I--il- -~--I- -~


-~~~~----
Li~-Tr.
aL--1_11~=-=1
----r--- --~-
-~- -






-~-~'
.--.r-.--


-----------i---


-;--, -


\0





bf

rIL.
..... C






77


:Ifsa -----/--
...-c~t~-



--=I-~


:-~ ::::
-c-~~--


-- ------i~_
~

-=-~


g
--2;L
9.
~-"














APPENDICES













APPENDIX A
JACOBI'S THEOREM


In evaluating hij in Equation (2.17) we used Jacobi's

theorem which is stated as

(A-1)(k)=IAj-ladj(k)A. (7.1)

The derivation of the theorem is quite lengthy, and so we

will add here only a brief explanation on the notions in

Equation (7.1) and the application of the theorem in evalu-

ating h... A reader interested in the development of the
theorem could refer to Aitken(2).
Let A be a matrix of order nxn. The determinant obtained

by suppressing m rows and m columns of IAlis called a minor

of IAl of order n-m.
(R)
The kth compound A(k) of A is obtained by forming its

elements with minors of IAI of order k in the following ways
let all minors which come from the same group of k rows (or

columns) of A be placed in the same row (or column) of this
derived matrix, and let the priority of elements in rows or
columns of this matrix be decided on the principle by which
words are ordered in a dictionary or lexicon. For example,
minors from rows 1,2,4 of A will appear in earlier rows than
those from 1,2,5 or 1,3,4 or 2,3,4; and similarly for columns.

Consider a minor IBI taken from rows il,i,2 im and








columns jl,j2, ...jm of A. Then the complementary minor

formed by the remaining n-m rows and n-m columns of A, with
the sign factor (-1)i'+'**+im+jl+'''+jm multiplied in front,

is called a cofactor of minor IBI.
Then the kth adjugate compound of A which is denoted by

adj(k)A is obtained in the following ways take the kth com-
(k)
pound A (k), replace every element in it by its cofactor in
IAI, and transpose the resulting matrix.
Now let the matrix M in Equation (2.8) be equal to A-1

in Equation (7.1). Then


(7.2)


and the determinant


(7.3)


is the element in the 1st row and 1st column of the nth

compound matrix of M. This element, if divided by IMI, should
be by the Jacobi theorem equal to the element in the 1st

row and 1st column of the nth adjugate compound of M-1, which

is


D2= mn+l,n+l...*mn+l,2n



m2n,n+1' m2n,2n


where the matrix (mij)=M-1. Thus
1;'


(7.4)


M(k)=IM-11 -ladj(k)M-1,







D1=IMID2=DD2. (7.5)
Since

I(hij) =D2-, (7.6)

by substituting Equation (7.5) into Equation (7.6) we obtain
the relation in Equation (2.16)

I(h ij =D/D1.

Now let


P'(Pij)= mn+1,n+1' ..,mn+1,2nN
(7.7)
m2nn+l '...... m2n,2n


Then, from Equation (2.11), the (i,j)th element of (hij) is

hij=D2-11Pjil, (7.8)

where JPjil is the cofactor of the element Pji, which is
obtained by deleting the jth row and ith column in IPI and
multiplying it by the sign factor (-1)i+
But this cofactor could be also obtained by suppressing
11
n+1 rows and n+1 columns from the 2nx2n matrix M-1. That is,
we suppress from M- the first m rows and m columns, the
(n+j)th row and (n+i)th column, and multiply the sign factor
(-1)i+j. By the Jacobi theorem, however, we can evaluate
this cofactor by evaluating the corresponding minor in M.
From Equation (7.2) we note that the corresponding minor
consists of the first m rows and m columns, the (n+i)th row








and (n+j)th column in M.
Thus


SPjil =l-1


R11l .... Rpn' Rl


R ,. R R '
nl'""' -in' Rnj
-Rill "*,* in '-R ij


and, substituting the above equation into (7.8), we finally
obtain as in Equation (2.17)


h ij=D2-1PjiI

=D 1 RLI." ., Rin' RIj'

n'
SRnl' '' R nn, R'
I i
_i R in ,i "


(7..10)


(7.9)













APPENDIX B
EVALUATION OF A DETERMINANT


Since the OS-8 operating system for the PDP-8/e computer

did not have in its library a program to compute a determi-

nant of a matrix, a subroutine is provided to invert a matrix

or evaluate a determinant of a matrix. The Gauss-Jordan

pivotal elimination method is employed here.

In the evaluation of a determinant by the pivotal elimi-

nation method, two rules of transformation for determinants

are involved

1). If a matrix C equals- the product of matrices A and

B, then ICI=IAIIBI. (7.11)

2). If a matrix B is formed-by interchanging row i and

row k of A, then JBI=-jAI. (7.12)

If a non-singular matrix A can be reduced to the

identity matrix I by premultiplying A by a matrix B

BA=I, (7.13)

then, by postmultiplying both sides by A-1
BA-1=IA-1

or

B=A-1, (7.14)

that is, we see that B is the inverse matrix of A.

The operation in Equation (7.13) could be carried out

in n transformation steps








E En-_...E2E1A (7.15)

successively diagonalyzing one column at-each step. For

example, at the first stage the resultant matrix is


1/a11,0,......... ,0

-a21/a11,1,0,.....,0

-a31/a11,,,0,,0...,0


11' a2 ..., ain

n'

anl,an2,...,ann
nj


1, a12/all,..., a n/all

O,a22-a21ai2/all, ...,an-a2lan/all



0,an2 -anl 12/all ,...,a-anlalr/a 1


=A(1)


(a11#0)


S (7.16)


In general, after k steps, we want

EkEk- ...E1A


k
1,.......,0,alk+1,..

0,1,0,...,0,a k+l,.
2, k+1'


k
',al,n

ia2,n
k.<
.

Sk k
0'00'',...,lak,k+1...' .k,n

k k
0,0,0,...,0,a k+l ...,a ,n


=A(k)


(7.17)


where the superscript k indicates that the elements belong

to A(k). Then, after n steps


(7.18)


E A=


.


EnEn-l '...EA=I.





85
(k)
In order to produce A(k the matrix Ek should have the
form
Ek= 1,0,0,...,-alk /akk ..., k-1
k-= k-1 (akk ~ 0)
0,1,0,...,-a2k /akk ,...,


0,0,0,..., 1/akk ,... (7.19)


k-1/a k-1
0,0,0,...,-a.k /a kk '',.


We note
|Ek =l/ak1. (7.20)

The operation in Equation (7.18) could be expressed with
recursion formulas which calculate the elements of A(k) from
the elements of A(k-i) .

ak k-1/ak-1 (a k-l,,
akj-akj /akk (kk 0)
a (j=1,....n), (7.21)
k k-1 k-1 k W
i=i -aik kj, (i
0
for k=l,...,n, provided we define a..=a...
The procedure to obtain the inverse matrix is then

EEnEn-...E1I=A-1. (7.22)

Let us denote

B(k)=EkEk-1. .E1 '
(k)
then the elements of B(k) are calculated from the elements
of B(k-1) by the recursion formulas








bk =bk-1/ak-1
kj kj "kk
S(j=1,...,n) (7.23)
bk =bk-1 ak-bk ,(ik)
1j= i -ik k\ '

Since, from Equation (7.18),

1=II =IEnEn-1...EIAI,

using the rule in Equation (7.11) and the result in Equation

(7.20), we obtain

1=IEn JE ln1 ... llj IA

=fan-1an-2 2 1 a --
n(a nan-l,n-l*a 33a22a11) AI'
or

Al=aj an1 (7.24)
AJ= 11 22 ...nn


When the computations are carried out by a digital
computer, the roundoff error occurs in each stage of the
elimination process and propagates to the next stage. An
analysis of the error expressions corresponding to the
recursion formulas for the elimination process reveals that
this roundoff error can be minimized by proper choice of the
ratio ail/akk (McCalla). To reduce the error, the numerically
largest element is selected as the element akk, which is
called a pivot element, from the elements aij (i,j ,k). In
order to bring the largest element to the pivotal position,
row and column exchanges are necessary and the corresponding
sign change in the determinant should be applied by the rule
in Equation (7.12).














APPENDIX C
NUMERICAL INTEGRATION


1). The error function in Equation (3.13)

erf(x)=2T2 fdtexp(-t 2) (7.25)
0


is evaluated by the rational approximation (Abramowitz)

erf(x)=l-(alt+a2t2+a3t3+a4t4+a5t5)exp(-x2)+r(x), (7.26)

where (0 t=(1+px)-1,
p=0.3275911,
a1=0.2548296,
a2=-0.2844967,
a3=1.4214137,
a4=-1.453152,
a5=1.061405.

The error bound in this expansion is

Ir(x)l( 1.5xio-7.

2). To evaluate the statistical function

Q(x)=(2v)-' Sdtexp(--t2), (7.27)
x

the following expansion is used (Abramowitz)s








Q(x)=Z(x)(bt+b2t2+b t3+b t4+b t5)+r(x), (7.28)

where (0
Z(x)=(2Tr) exp(-x2 ),
1
t=(l+px)1,
p=0.2316419,
b1=0.3193815,
b2=-0.3565638,
b =1.7814779,
b4=-1.821256,
b =1.330274.
5-

The error bound is

Ir(x)l| 7.5x108.

For a large value of x, the approximation of Q(x) by
continued fractions

Q(x)=Z(x)[1/(x+l/(x+2/(x+3/(x+...)))], (x>0) (7.29)

would be a better evaluation. But, since the difference
between the two approximations was not significant enoughto show
up on the plotted graphs, the expansion in Equation (7.28) is
used in the computation throughout the range of x.

3). For other numerical integration, the trapezoidal
rule is employed.
Suppose f(x) is defined at equally spaced points xi
(i=0,1,...,n), then the area under the curve of f(x) between

x0 and xn is computed by the trapezoidal rule
n-1
T= E Ih[f(xi)+f(xi+h)], (7.30)
i=O






89

where h is the distance between two points of x

h=xi+1-xi, (i=0,1,...,n-1). (7.31)

Writing the above formula in expanded form and substi-
tuting the relation xi+l=xi+h, we obtain the classical form
of the trapezoidal rule

T=h[f(x2f()+2f(x)+2f(x)+...+2f(xn)+f(x )]. (7.32)













APPENDIX D
INTERPOLATION BY SPLINE FUNCTION


When an interpolation of data is needed, the spline

function is used to approximate the intermediate values of

the data. The spline function is a numerical representation
of the curve of a flexible strip or spline used in drafting

to draw a smooth line through a set of points (Carnahan).
Suppose we have n functional values (xi,f(x ))
(i=l,...,n), and want to interpolate the data over the

interval (x1,xn). The restraints that define the spline
function over interval (x1,x ) are then as follows:
1). In each interval the spline function is a cubic poly-
nomial, that is

P3,i(x), (i=1,2,...,n-1) (7.33)

where the subscript in p3 indicate that p is a polynomial of
degree 3.

2). Each cubic polynomial p3,i(x) should match the functional
values at its two ends x. and xi+l


P ,i(xi)=f(xi)'
S (i=i,2,...,n-1) (7.34)
P3,i (xi+)=f(xi+),

3). The first and second derivatives of successive polynomials
must match each other in value at the intermediate points:





91

3,'i(xi)P i-'(xi)'
S(i=2,3,...,n-1), (7.35)
P3,i(xi)=P3,i-1(x ),

and the second derivatives vanish at two end points:

,1(x )=0,
(7.36)
3 ,n-1( n)=

Equations (7.34) through (7.36) amount to 4(n-1)
relations, which enables us to determine the 4(n-1) coeffi-
cients needed to completely describe the spline function in
Equation (7.33).
For convenience, define

f=f (x ),

h =x+1-xx
hi =Xi+l-xi' (7.37)

gi=P" 3,i-(xi)=P"3 i(xi)

Since each p" (x) is a linear function of x, it can be
3,1
obtained on the interval (Xi,xi+) by linear interpolation
of the values gi and gi+1 at each end of the interval:

,i(x)=gi(xi+1-x)/higi+1(x-xi)/hi' (7.38)

By integrating twice and imposing the conditions in Equation
(7.34), we find that

P3,i(x)=g (xi+l-x)3/(6hi )+gi+(x-x)3/(6h )+


(7.39)


(fi+A/hi-gi+lhi/6)(x-xi)+








(fi/hi-gihi/6)(xi+1-x), (i=1,2,...,n-1).

By differentiating Equation (7.39) and imposing the
conditions in Equation (7.35), we obtain the following system
of linear equations in the gii


gi-_lhi-i/hi+2gi(l+hi- /hi)+gi+l

=6[(fi+l-fi)/hi-(fi-fi-)/hi-i]/hi' (7.40)

(i=2,3, ...n-l),
and

gl=0,
gn=0, (7.41)

which follow from Equation (7.36).
Note that Equations (7.40) and (7.41) amount to n
simultaneous linear equations in the n unknowns gi1, ... n
Thus the values of gl,...,gn can be determined, and the
substitution of the obtained gi (i=1,2,...,n) into Equation
(7.39) yields the individual cubic spline polynomials for
use over successive intervals.














APPENDIX E
LIST OF PROGRAMS




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs