Full Citation 
Material Information 

Title: 
A signal detection from noise utilizing zerocrossing 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 ) nonfiction ( marcgt ) 
Notes 

Statement of Responsibility: 
by Woon Cheon Yeo. 

Thesis: 
ThesisUniversity of Florida. 

Bibliography: 
Bibliography: leaves 132135. 

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 

Full Text 
A SIGNAL DETECTION FROM NOISE
UTILIZING ZEROCROSSING 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. ZEROCROSSING INTERVAL DENSITY FUNCTIONS 8
Derivation of Zerocrossing Interval Density
Function . . . . . . . . . . 8
Evaluation of W for Gaussian Process . . . 12
Zerocrossing Interval Density Function for a
Bandlimited White Gaussian Noise . . . . 20
Evaluation of W for Sine Wave Plus Noise . . 23
Zerocrossing 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 ZEROCROSSING 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
zerocrossing 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
zerocrossing intervals of the process has been found to be
very efficient in detecting signals from noise. The operating
characteristic of the zerocrossing 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 zerocrossing measuring detector.
The operating characteristic of the zerocrossing
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 signaltonoise 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 zerocrossing 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 zerocrossing 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 zerocrossing measuring detector,
To compute the detection and false alarm probabilities of
the zerocrossing measuring detector, we need the zerocross
ing interval density functions of the processes with noise
alone and signal plus noise. The distribution functions for
the successive zerocrossing intervals of the bandlimited
Gaussian noise have been developed by Rice(36,37), McFadden
(26,27) and LonguetHiggins(23,24). Experimental works on
S.
the zerocrossing interval distribution of Gaussian noise
have been presented by Rainal(32) and Blotekjaer(6).
Rice's zerocrossing interval density function agrees well
with experimental data for smallvalues of zerocrossing
interval. The zerocrossing 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 LonguetHiggins fits closely to the experimental
data given by Blotekjaer. LonguetHiggins 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 LonguetHiggins approximation is used
here to compute the zerocrossing interval distribution of a
bandlimited Gaussian process.
In the signalplusnoise case, the zerocrossing inter
val density function is approximated by the first term of
the LonguetHiggins 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 zerocrossing 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 zerocrossing interval density for sine
wave plus noise in this study is based on the treatment of
Cobb.
Derivations of the zerocrossing 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 LonguetHiggins' 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
zerocrossing interval density for signal plus noise are also
presented in Chapter II.
In Chapter III, the schematic of the zerocrossing
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 zerocrossing
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 PDP8/e. Without an extended
arithmetic unit, the speed of the PDP8/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 zerocrossing interval density functions
and evaluated the statistical functions from their series
expansions.
CHAPTER II
ZEROCROSSING INTERVAL DENSITY FUNCTIONS
Derivation of Zerocrossing 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 upcrossing or downcrossing to indicate a crossing
with f'(t,)>0 or f'(t1)<0 respectively.
Consider the instance when a random function f(t) has
upcrossings 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 downcrossings 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 upcrossing in the interval (tl,t.+dtl),
downcrossing in (t2,t2+dt2), and an upcrossing 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 upcrossing at
t=t1 and a downcrossing 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 upcrossing at t=tl and a downcrossing after an interval
r, containing 2n crossings anywhere in between the two
crossings.
If f(t) has an upcrossing at t=t1 and a downcrossing
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 upcrossing at t=t1 and a downcrossing at t=t2,
is then the sum of the probabilities of all possible cases
when f(t) has an upcrossing at t=t1 and a downcrossing
after an interval T=t21t, 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 upcrossing at t=t1 and a downcrossing at t=t ,
containing a downcrossing at t=t2 anywhere in between t1 and
t Since the crossings at t=tI and t=t3 are an upcrossing
and a downcrossing respectively, f(t) can only have an even
number of crossings in the interval (tl,t3), and one of the
downcrossings among these crossings should pass zero at t=t2.
Suppose we have 2n crossings in (tl,t3) where t3tl=tr Then
the downcrossing at t=t2 could be one of those n down
crossings, and, therefore, there are n different possibilities
of f(t) having a downcrossing at t=t2. Since the downcross
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 upcrossing at t=t1
and a downcrossing 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 upcrossing at t=tI and a downcrossing at t=t3
containing at least one downcrossing 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
zerocrossing interval density (LonguetHiggins):
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 upcrossing 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 upcrossing 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
upcrossings 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(tti)=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+)]
=6R(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(tr) f(t)l
=Ev[6f (tjT)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)]
6r i
=E[ FAf(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(iXTMlX)
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=M1. 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)lexp(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 (LonguetHiggins)
W(+,...,+)=(2T)nD1 (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)PD12h11J 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) [(1v12) 1Vcos1O(12)], (2.21)
J3=(2n)3/2[ (vij)lI+(slbi+s2b2+s3 3)], (2.22)
where
sl=cos v31 v12v23)(1v312)(1v2 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 (R11R22R12R21) (h11h22) (1v122+
v12cos1 (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 downcrossings
in W. Suppose the kth zerocrossing is downcrossing. 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 M1, and hence the kth
row and column of (hij) and (vi ). Thus by changing sign of
v12. in Equation (2.23) we obtain (LonguetHiggins)
W(+,)
=(2)2(R11R22R12R21)(h1h22 v1212cos v12],
(2.25)
and also by changing corresponding signs of v.. in s and b
in Equation (2.23), we get
W(+,,)=(2n)3Dl (hh22h33) [I(v) +sb +
, 11 22 33 j~ U 1b
(2.26)
Zerocrossing 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)=flsinr.
(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 zerocrossing 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=t3t.l
The computation repeats with different value ofr 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 zerocrossing 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 signaltonoise ratio
is fairly high so that the zerocrossing intervals are mainly
determined by the signal frequency, and the probability of
having two downcrossings 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)=7lsinr,
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)=gis(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'=glsi')
=(2y)1Mr 'exp( XTM1X), (3.5)
where
X= s ] = acose
g1s1 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 expCiCa2cos2 +
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)=gis'(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)
glS1 gl+qasine
g2s21 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
Lg2+qasin(8+qq )
The exponent in Equation (3.10) may be written as
M' 1
iU'M U
=(21MI) 12a2M11[cos2e+cos2(e+qr)]+2a2M14cosecos(e+qT)
2aM12[ (gl+qasine )cose+[g2qasin(e+qr)]cos(e+qr)]
+2aM13 [gqasin(o+q)]cose+(glqasin )cos(++qsr)]
+M22[ (91+qasine)2+[g2qasin.(e +q)]2]
2M2 (g1+qasine) g2qasin(9+qr)] (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")(1R2)(R')2
M23=(R")(1R2)+R(R')2,
!MI =[ (1R)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 =M22M23cos(2x)
=(M22M23)cos2 +(M22+M23 )sin2 X
D2=alMI ([[ (M12M13 )cos(qr)+q(M22M23) sin (iq)])
cos(e+Fqr)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(M12M13)sin(iqT)cos(GqT)
+2 (M22M23) sin2 r) ]cos2 (+1qr)
+[(M11M14)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,=1R,
C2=1+R,
C =RO"+R",
29
C4=R0 "R",
C=(IR)(RR")(R')2=01C 4002,
C6=(1+R)(R"+R")(R' )2 =C2CC2,
=C7COcs(iqT)+qC2sin(iqr),
C8=C0sin(iqT)+qC1'cos(iqT),
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]
+C11 (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 dxD12cos(2X)[(l+D 2)exp(D3)
I r 2
nD4 (3/2+D42)(1erfD )exp(D5 ),
(3.13)
where
D4=D2/D1 ,
D=D 3D 2
Zerocrossing Interval Density Function for
Sine Wave Plus Noise
The zerocrossing 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 38. 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 signalplusnoise 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 downcrossing between an upcrossing and a downcrossing
separated by an interval T=t t,. When the amplitude of the
signal is high, the zerocrossing intervals of the signal
plus noise are close to the half period of the signal wave,
Th=n/q, and the probability that the zerocrossing 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 zerocrossing 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 zerocrossing 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 zerocrossing
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 zerocrossing 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
....... ..... .. ....... ...... ......... . . .
Ug
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
signaltonoise ratio and the length of a signal are known
to us. We also assume that the signaltonoise ratio is fairly
large so that the additive noise on a signal only disturbs
the zerocrossing interval determined by the signal frequency
and does not add more zerocrossings in the process. Then we
know the number of zerocrossings needed to observe the signal
in full length.
The schematic of the detector is shown in Figure 9.
Zerocrossipgs of f(t), both upcrossings and downcrossings,
are first detected and the output pulse of the zerocrossing
detector is used to strobe the system. The intervals between
zerocrossings 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 zerocrossing
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 zerocrossing
interval comparator.
The outputs of the zerocrossing interval comparator are
stored in a shift register and the content of the shift regis
ter is counted by an updown counter. If the signal s(t) is
composed of n cycles of sinusoidal waves, then there are 2n
zerocrossings 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 i2n' (4.2)
ci= (4.2)
1 Ci +1, if Z =1 and Zi2n=0
Ci_1l, 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=i2n
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=0H1)=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 2nk (5)
k=m s s
The corresponding false alarm rate is
Pf=P(D =1 H )=P(C imH )
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~ure 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 NeymanPearson 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 )2nk (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(lt)nm]/[ dttm1(lt)nm (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 zerocrossing interval density
function of noise from T=0 until the integrated area equals
Pn, and the zerocrossing 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
zerocrossing 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 1519. 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=102 is 6, while 7 is the
optimum for Pf=10 .
45
U ..
..... ............................... ...
.O
0o 0 r
a1
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 n1 ~
33
. ... ...... ..... ........ .......... ... ........ ....... .. ... .... ...
........ ........ ...... ...... . ... ..... ..... .
............. . ... .. ........ ....... ........ ........ ........ ...... ...... .
.......... ........ ........ ........ ........ ....... ........ ...... ....... ....... ..
:o o o o o o
'.4
*r4
ae
....*. . . . .. ......... . ............ ....... .... ... .
............ ........ ................. ........ .......... ........ ....... ....... 0
................ . ........ 1 ........ ........ ........ ...... ........ ... .... ..
'Cr IrTIii i1ir*
 o m I I I?
o CCoo N cooo 
4.
I :
'.4
I . ~.......,.......i....... ..
1  r rrr i rr
0 S fl Y f .
.5 a a C a 0 0
acr o n
46
c
a
......,..... ....... ....... ................. .
... ............... ........ .. ..... ........ ....... ....... ........
o
........ ............... ........ ....... ....... ........ ....... .
0
.. ..... ..... .. ............. . .................
................. ................ ....... ......... ...... .. ........ .....
............... ......... ......... ................ .......... ................. ........
o:0
C(

0
I Cf
"""""""~""""""""" ""~" ~a
 0 0 0 0 0 0 0 0 0
OcrecDrwnnr.
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 rr 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
fla
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
00C *
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(Tt) 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 NeymanPearson 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 TM1X), (5.6)
where X is the column vector of sampled values. By writing
Q=M1.(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.kshk (h=0,l), the above equation could be
written as
p(n)=p(fsh)
=(2T)mlMJj expL[ q(f sihi jfshj
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 (fksk)
(k=l,...,m) is
PO(f)=p(fs0)=p(n)
1 1
=(2n)MIMexp[ E qij(fioi)(fjoj)].
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(fsl)=p(n)
1 1 2
=(2)_ImlJMI exp[i qi isli )jf 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)(fslj)]/
m
exp[i E qi.(f soi)(f s
i, =1 J a
m
=exp[j E qj(2f s0j2f sljsOi1s+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 jt)= 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[ (kslk)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 ksk)2/exp[i 2 (fSok) J (5.17)
k=1 k=
or, upon combining terms,
r(f)=exp[2 E (2fks0k2fks ks0k +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(kslkf ksOk+sk ~slk )n(ro). (5.19)
k=l1
Let
m
G= E (f kskfksOk Ok2s 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+S0kji lk )]
k=l
Since Sok=O,
m
Eo[G]=1 8sk2
k=1
The variance of G under H0 is
Vo[G]=E[(GEo[G])2].
But
GEO[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) Sexp[(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
GE1[G]= E (slk+nk)slk+ (s ) 12
k=1 k=l k=1
m
= E nksik.
k=1
VI[G]=E[(GE1[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(GiS)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 (GS) 2/]
c
1 0 1 2
=(2n)S ~J' dzS exp(iz2)
(cIS)/S2
=(24) dzexp(iz2) (5.31)
1
(cIS)/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
zerocrossing 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 4x107 while the false alarm
probability of the zerocrossing measuring detector at the
same detection rate is 3x103.
This superior performance of the correlation detector
to the zerocrossing 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"  ,TT ',Fj
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 zerocrossing 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 zerocrossing measuring detector.
S66
*b * 1 *
. .
_______________ S 
45
_ _ _ _ _  ^
 4 ?:ySsr
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[(GE2[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)SexpC(GE2) /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(GE2 2/S]
c
=(2T) X S dzexp(iz2), (6.5)
1
(cE2)/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 zerocrossing
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
zerocrossing 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 zerocrossing 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 zerocrossing 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 zerocrossing 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 zerocrossing 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.
aa
.....1 :......... ....l... ........ ......... ........ ..
L,0
C0
...... ......
....... .... .. .. ... .... ........ ....... ..... .. ... ... ..... ........ j
00 0 u0 0
*. .* '.. . .. .. .
.. ....... .......... ....... .' ........ ....... I ......... .
3 0 ,
...... ........ .... ........ .. .... ... ....... .... . .......... . .i. . .
1 I
. .. . .. . . . . .
 "
as s1 and the signaltonoise 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 zerocrossing 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 zerocrossing 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[(GE[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)4Ih1S .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 zerocrossing
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 zerocrossing 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
zerocrossing 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 zerocrossing 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
zerocrossing 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 zerocrossing 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~~
"7ie
Wc~
~a
=~;5=
4~t
~CL ~
~:~S~L
L  T;
i
~~. .~~r~;r_
`i: ~c
L_~L~
 ~~;~~
,

   r~s
 ~~1 ____
Iil ~I ~
~~~~
Li~Tr.
aL1_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
(A1)(k)=IAjladj(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 nm.
(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 nm rows and nm 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 A1
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 M1, which
is
D2= mn+l,n+l...*mn+l,2n
m2n,n+1' m2n,2n
where the matrix (mij)=M1. Thus
1;'
(7.4)
M(k)=IM11 ladj(k)M1,
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=D211Pjil, (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 M1. 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 =l1
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=D21PjiI
=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 OS8 operating system for the PDP8/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 GaussJordan
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 formedby interchanging row i and
row k of A, then JBI=jAI. (7.12)
If a nonsingular 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 A1
BA1=IA1
or
B=A1, (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 ateach 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,a22a21ai2/all, ...,ana2lan/all
0,an2 anl 12/all ,...,aanlalr/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=
.
EnEnl '...EA=I.
85
(k)
In order to produce A(k the matrix Ek should have the
form
Ek= 1,0,0,...,alk /akk ..., k1
k= k1 (akk ~ 0)
0,1,0,...,a2k /akk ,...,
0,0,0,..., 1/akk ,... (7.19)
k1/a k1
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(ki) .
ak k1/ak1 (a kl,,
akjakj /akk (kk 0)
a (j=1,....n), (7.21)
k k1 k1 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=A1. (7.22)
Let us denote
B(k)=EkEk1. .E1 '
(k)
then the elements of B(k) are calculated from the elements
of B(k1) by the recursion formulas
bk =bk1/ak1
kj kj "kk
S(j=1,...,n) (7.23)
bk =bk1 akbk ,(ik)
1j= i ik k\ '
Since, from Equation (7.18),
1=II =IEnEn1...EIAI,
using the rule in Equation (7.11) and the result in Equation
(7.20), we obtain
1=IEn JE ln1 ... llj IA
=fan1an2 2 1 a 
n(a nanl,nl*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.5xio7.
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
n1
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+1xi, (i=0,1,...,n1). (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,...,n1) (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,...,n1) (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,...,n1), (7.35)
P3,i(xi)=P3,i1(x ),
and the second derivatives vanish at two end points:
,1(x )=0,
(7.36)
3 ,n1( n)=
Equations (7.34) through (7.36) amount to 4(n1)
relations, which enables us to determine the 4(n1) coeffi
cients needed to completely describe the spline function in
Equation (7.33).
For convenience, define
f=f (x ),
h =x+1xx
hi =Xi+lxi' (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+1x)/higi+1(xxi)/hi' (7.38)
By integrating twice and imposing the conditions in Equation
(7.34), we find that
P3,i(x)=g (xi+lx)3/(6hi )+gi+(xx)3/(6h )+
(7.39)
(fi+A/higi+lhi/6)(xxi)+
(fi/higihi/6)(xi+1x), (i=1,2,...,n1).
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_lhii/hi+2gi(l+hi /hi)+gi+l
=6[(fi+lfi)/hi(fifi)/hii]/hi' (7.40)
(i=2,3, ...nl),
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

