REALTIME COMPOSITE SIGNAL DECOMPOSITION
By
DAVID PRESTON SKINNER
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
Dedicated to Betty became I tove heA.
__
ACKNOWLEDGEMENTS
The author wishes to thank Dr. D.G. Childers, chairman of his
supervisory committee, for his guidance during the preparation of
this dissertation. The example set by him, both intellectually and
personally, has been a highlight of this author's experience at the
University of Florida. The author also wishes to express his gratitude
to the other members of his supervisory committee for their assistance.
Finally, the author would like to express his appreciation to his
parents, Mr. and Mrs. J.P. Skinner, for their encouragement and
love throughout the years, and to his wife for her unflinching support.
____
TABLE OF CONTENTS
PAGE
ACKNOWLEDGEMENTS ......................................... iii
LIST OF TABLES ....................................... .......... vii
LIST OF FIGURES ................................... ..... .. viii
ABSTRACT .............................................. .... xii
CHAPTER I INTRODUCTION ................................... 1
CHAPTER II THE CEPSTRA .................................... 5
THE POWER CEPSTRUM ........................... 5
History ....................................... 5
Introduction ............................... 6
THE.COMPLEX CEPSTRUM AND WAVELET RECOVERY .... 10
History .................................... 10
Introduction ............................... 11
The Relation Between the Complex and Power
Cepstrum ................................... 18
The Phase Cepstrum ......................... 20
SUMMARY ..................................... 21
CHAPTER III THE PITFALLS OF CEPSTRUM COMPUTATION ........... 23
LINEAR PHASE TERMS ........................... 23
PHASE UNWRAPPING ERRORS ...................... 43
ALIASING ..................................... 47
OVERSAMPLING ................................. 49
SUMMARY ...................................... 49
(iv)
PAGE
CHAPTER IV THE EFFECTS OF SOME DATA PROCESSING TECHNIQUES
ON THE CEPSTRA ................................. 51
WINDOWING THE COMPOSITE SIGNAL ................ 51
The Exponential Window ..................... 53
Interpretation of Results .................. 55
The Common Window .......................... 58
Interpretation of Results .................. 61
Conclusions ................................ 76
The exponential window ............ ....... 76
The common windows ... ..................... 77
WINDOWING THE LOGSPECTRUM ................... 77
Interpretation of the Results .............. 80
Conclusions ................................ 81
HANNING THE LOG SPECTRUM ..................... 81
Experimental Results ....................... 84
Conclusions ............................... 86
THE ADDITION OF ZEROES ....................... 87
Experimental Results ....................... 91
Conclusions ................................ 110
CHAPTER V ALGORITHMS FOR A REALTIME WAVELET RECOVERY
SYSTEM ........................................ 112
THE DFT ALGORITHMS ........................ ... 113
The CooleyTukey (decimation in time) and
SandeTukey (decimation in frequency)
Algorithms..... ............................ 113
Modifying the Basic FFT Algorithms ......... 114
The Bergland Algorithm ..................... 118
(v)
PAGE
The Hartwell Modification .................. 118
Which Algorithm for Complex Cepstrum
Computation? ................................ 122
COMPUTATION OF NONLINEAR FUNCTIONS ............. 126
Computation of Nonlinear Functions .......... 126
PHASE UNWRAPPING AND LINEAR PHASE REMOVAL ..... 131
LINEAR FILTERING ... .......................... 132
:AN OVERALL LOOK AT SYSTEM PERFORMANCE ......... 133
An.Estimate of System Performance ........... 136
REALTIME COMPUTATION OF THE POWER CEPSTRUM ... 136
SUMMARY ...... ................................. 137
CHAPTER VI SUMMARY AND CONCLUSIONS ......................... 139
COMPUTATIONAL PROBLEMS ........................ 139
WINDOWING ..................................... 140
EXTENDING THE DATA RECORD WITH ZEROES ......... 142
ALGORITHMS FOR REALTIME WAVELET RECOVERY ..... 142
SUGGESTIONS FOR FUTURE RESEARCH ............... 143
APPENDIX A A COMPARISON OF THE ECHO DETECTION CAPABILITY
OF THE PHASE, POWER AND COMPLEX CEPSTRA ......... 145
APPENDIX B THE EFFECTS OF ADDITIVE NOISE ON THE CEPSTRA
AT HIGH SNR ................................. .... 155
APPENDIX C THE INVERSE TRANSFORM OF THE TRANSFORM OF A
REAL VALUED SERIES UTILIZING THE FORWARD
ALGORITHM .................................. ... 173
APPENDIX D PROGRAM LISTING ................................. 174
BIBLIOGRAPHY ................................................. 185
BIOGRAPHICAL SKETCH ............................. ........... 187
(vi)
LIST OF TABLES
TABLE PAGE
1 Algorithms for Cepstrum Computation ................ 123
2 FFT Processor Performance ............................ 125
3 Throughput Comparison of Stages in the Wavelet
Recovery System ................................... 135
(vii)
LIST OF FIGURES
FIGURE PAGE
1 Computation of the Power and Complex Cepstrum ....... 8
2 Phase Unwrapping .................................. 14
3 Effect of a Delay in the Composite Signal on the
MSE of the Recovered Wavelet ...................... 25
4 Composite Signal ................................... 28
5 Unwrapped Phase Curve .............................. 29
6 Log Magnitude ...................................... 30
7 Complex Cepstrum ................ ......... ........ 31
8 Phase Cepstrum ..................................... 32
9 Power Cepstrum ................... ................. 33
10 Recovered Wavelet .................................. 34
11 Composite Signal .................................... 36
12 Unwrapped Phase Curve ... ..................... ... .37
13 Log Magnitude ....................................... 38
14 Complex Cepstrum ............................... 39
15 Phase Cepstrum ....................................... 40
16 Power Cepstrum ....................................... 41
17 Recovered Wavelet .................................... 42
18 Phase of X(ejWT) .................................... 45
19 Phase Unwrapping Errors ........................... 45
20 MSE of Recovered Wavelet when the Input Data Record
is Exponentially Windowed ......................... 55
(viii)
FIGURE PAGE
21 MSE of the Recovered Wavelet when the Input Data
AreHamming Windowed ................................. 60
22 Composite Signal ................................... 65
23 Unwrapped Phase Curve .............................. 66
24 Log Magnitude .................................. .... 67
25 Complex Cepstrum ................................... 68
26 Phase Cepstrum ..................................... 69
27 Power Cepstrum .................................... 70
28 Recovered Wavelet ................................... 71
29 MSE of the Recovered Wavelet when the Log Spectrum
Is Hamming Windowed ................................ 79
30 MSE of Recovered Wavelet when the Log Spectrum Is
Hanning Smoothed ................................... 85
31 Composite Signal ................................. . 93
32 Unwrapped Phase Curve ............................... 94
33 Log Magnitude ....................................... 95
34 Complex Cepstrum .................................... 96
35 Phase Cepstrum ..................................... 97
36 Power Cepstrum ...................................... 98
37 Recovered Wavelet ................................... 99
38 Composite Signal (Zeroes Added, 512 Points) ......... 101
39 Unwrapped Phase Curve ............................... 102
40 Log Magnitude ....................................... 103
41 Complex Cepstrum .................................... 104
42 Phase Cepstrum ...................................... 105
43 Power Cepstrum ...................................... 106
44 Recovered Wavelet .................................... 107
(ix)
FIGURE PAGE
45 MSE of the Recovered Wavelet when Data Record Is
Extended with Zeroes ................................ 109
46 CooleyTukey FFT Algorithm ......................... 115
47 The SandeTukey FFT Algorithm (with bit input
reversed) ......................................... 116
48 Bergland Real Valued Input Algorithm ................ 119
49 Bergland Inverse Algorithm ............... ......... 120
50 Computation of log(x) .............................. 128
51 Computation of exp(x) ............................. 129
52 Overall Wavelet Recovery Algorithm .................. 138
53 Composite Signal .................................. 148
54 Unwrapped Phase Curve .......................... 149
55 Log Magnitude ..................................... 150
56 Complex Cepstrum .................................. 151
57 Phase Cepstrum ..................................... 152
58 Power Cepstrum ............. ........................ 153
59 Recovered Wavelet .................................... 154
60 Composite Signal ................................. 158
61 Unwrapped Phase Curve .............................. 159
62 Log Magnitude ..................................... .. 160
63 Complex Cepstrum ................................... 161
64 Phase Cepstrum ..................................... 162
65 Power Cepstrum ......................... ........... 163
66 Recovered Wavelet .................................... 164
.67 Composite Signal (Zeroes Added, 1024 Points) ........ 166
68 Unwrapped Phase Curve .............................. 167
(x)
FIGURE PAGE
69 Log Magnitude ....................................... 168
70 Complex Cepstrum .................................... 169
71 Phase Cepstrum ..................................... 170
72 Power Cepstrum ..................................... 171
73 Recovered Wavelet ................................. 172
(xi)
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
REALTIME COMPOSITE SIGNAL DECOMPOSITION
By
David Preston Skinner
December, 1974
Chairman: Donald G. Childers
Major Department: Electrical Engineering
The purpose of this research is to investigate the recovery in
realtime of an unknown wavelet from a composite signal in the presence
of additive noise. Specifically, the composite signal consists of
a wavelet and its echoes.
Two techniques, the power and complex cepstra, are used in the
wavelet recovery procedure. Historically, these techniques have been
computed separately. An alternate definition of the power cepstrum
is presented which closely relates the two techniques,and makes
computation of the power cepstrum from the complex cepstrum trivial.
This unification led to the discovery of a third cepstrum technique,
the phase cepstrum, which like the power cepstrum is useful in the
determination of echo epoch times.
Several problems inherent in cepstrum computation which must be
overcome to achieve satisfactory wavelet recovery are examined. These
i.) 
problems include linear phase terms, phase unwrapping errors, aliasing,
and oversampling.
The effects of windowing (as used to reduce leakage in spectral
analysis) at various stages in the wavelet recovery system are examined.
In general, windowing of the input data or log spectrum is detrimental
to wavelet recovery. A windowing of the complex cepstrum may result
in some improvement in wavelet recovery.
The addition of zeroes to the input data record is explored.
This is observed to reduce the phase unwrapping and aliasing problems.
Finally, algorithms suitable for use in a realtime wavelet
recovery system are examined. A system based on the selected algorithms
should be able to achieve sampling rates of around 10 samples/sec.
(;;iii)
CHAPTER I
INTRODUCTION
Composite signals consisting of the convolution of two or more
simple signals arise in a number of physical situations. The
separation of the contributions of the components of such signals
often yields important information about the underlying physical
processes. Several techniques have been devised for this purpose,
these include inverse filtering [1], Wiener filtering [1], decision
theory [2] and the power and complex cepstrum techniques to be
discussed herein. Each of these techniques is applicable under its
own circumstances.
In particular the power and complex cepstrum have been shown
to be invaluable when the composite signal consists of a basic
wavelet convolved with a train of impulses [3]. One example of this
is a composite signal consisting of a basic wavelet plus echoes
as might arise in sonar or radar applications. The power cepstrum,
defined as the power spectrum of the logarithm of the power spectrum
of the given signal, has been shown to be successful in the deter
mination of echo arrival times when the composite signal consists of
a basic wavelet and its echoes [3]. Of course, since the power
cepstrum is derived from the power spectrum there is no hope of
reconstructing the original wavelet from it. The complex cepstrum
technique, defined as the Fourier transform of the complex logarithm
of the Fourier transform of the given composite signal, overcomes
this limitation since the phase information is preserved. The success
of both techniques stems from the fact that convolution in the
time domain is equivalent to multiplication in the frequency domain.
Thus by taking the logarithm of the frequency domain function the
contributions from the basic wavelet and impulse train may be
separated into a sum. In the case of the complex cepstrum appropriate
filtering techniques may then achieve separation of the components,
and by performing the inverse operation (to the complex cepstrum)
we may extract the basic wavelet. There is, in fact, no other
technique either linear or nonlinear, which allows the extraction
of an unknown wavelet from a composite signal consisting of the
wavelet and its echoes.
This investigation will primarily be concerned with the cepstrum
techniques as applied to echo detection and extraction. However,
the results are in most cases equally applicable to the other areas
in which cepstrum techniques have found applicability. Among these
areas are seismology [4,5], and speech [6,7].
In the past the cepstrum techniques have been limited to non
realtime applications. It seems clear that the performance of
these techniques in realtime would aid the research effort in the
abovementioned fields, and perhaps open new areas of interest.
This investigation is the first into the area of realtime compu
tation of the power and complex cepstrum, and wavelet recovery. This
is a formidable problem since these techniques involve the.computa
tion of several Fourier transforms together with nonlinear operations.
The purpose of this research was:
(1) to identify and investigate the problems inherent in
computation of the cepstra, and in wavelet recovery,
(2) to investigate the use of ordinary data handling techniques
(windowing and the addition of zeroes) at various points in
the wavelet recovery algorithm to improve the performance of
these techniques both in noise free, and additive noise
environments, and
(3) to select algorithms suitable for the realtime computation
of the cepstra, and for realtime wavelet recovery, and finally,
to estimate the possible data rate of a system based on the
selected algorithms.
Chapter II reviews the power and. complex cepstrum techniques.
Historically these two techniques have been treated separately. An
alternate definition of the power cepstrum is presented which closely
relates the two techniques, and makes computation of the power
cepstrum from the complex cepstrum trivial. This unification of the
two techniques led to the discovery of a third cepstrum technique,
the phase cepstrum, which like the power cepstrum is useful for the
determination of echo epoch times. This technique is introduced and
discussed. A comparison of this new technique with the complex and
power cepstra for the estimation of the echo amplitudes, and epochs
in the presence of noise is reported in Appendix A.
Chapter III presents some problems of cepstrum computation which
must be overcome to achieve satisfactory wavelet recovery. Some of
these have been reported by other authors ; others were previously
unmentioned. Problems of linear phase terms, phase unwrapping errors,
aliasing, and oversampling are discussed, and methods for their
alleviation are indicated.
Chapter IV analyzes (both theoretically and experimentally)
the effects of windowing at various stages in the wavelet recovery
process, and the effects of extending the input data record with
zeroes. Considerable insight is gained into the performance of
cepstrum techniques in additive noise.
Chapter V discusses the algorithms suitable for use in a
realtime wavelet recovery system, and roughly estimates the per
formance of such a system. The computation of the power cepstrum
(only) in realtime is also discussed. And, finally, Chapter VI
collects and summarizes the results of the previous chapters.
CHAPTER II
THE CEPSTRA
This chapter presents a brief history of, and introduction to,
the power and complex cepstrum techniques. In addition, a new
cepstrum technique, the phase cepstrum, is introduced along with
an alternate definition of the power cepstrum which serves to unify
the cepstra, and is of computational interest.
THE POWER CEPSTRUM
History
The power cepstrum was first presented by Bogert et al..[4] in
1963 as a heuristic technique for finding the echo epoch times of
a composite signal. Subsequently, A.M. Noll [6] in 1964 successfully
applied the power cepstrum to speech data to determine the pitch
period and for voicedunvoiced detection. Noll also simulated a
technique for short time power cepstrum computation utilizing the
direct computation of the discrete Fourier transform (DFT). In 1966,
Bogert and Ossanna[8] examined the statistical properties of the
power cepstrum of a Gaussian signal in Gaussian noise. Halpeny [9]
examined the properties of the cepstrum of a composite signal
consisting of a basic wavelet plus multiple echoes in the presence
of additive noise. Kenerait [3] then extended Halpeny's results
and examined the effects of echo truncation on the power cepstrum
(in conjunction with his work on the complex cepstrum). Kemerait's
dissertation shows the power cepstrum to be invaluable in the
detection and estimation of echo amplitudes and epochs.
Introduction
Basically, the power cepstrum, usually referred to in the
literature as simply the cepstrum, is just a clever way of separating
the contributions of two simple signals to the power spectrum of
a composite signal (which is the convolution of the two simple
signals). Of course to effectively separate any two signals there
must be some difference on which we may base the separation; for
the power cepstrum (and analogously for the complex cepstrum) this
difference is that the power spectra (more properly the logarithms
of the power spectra) of the two simple signals must vary (as a
function of w) at different rates; i.e., in Tukey's terminology [4]
they must occupy different ranges of quefrency. Historically the
power cepstrum has been defined as the power spectrum of the logarithm
of the power spectrum of the given signal. In actuality, the power
spectrum of the logarithm of the power spectrum does not exist for
most signals. The power cepstrum is only meaningful when defined
in a sampled data sense (as is the complex cepstrum). Thus the
following definition is proposed: The power cepstrum is the square
of the inverse ztransfonn of the logarithm of the magnitude squared
of the ztransform of the given sequence. Evaluated on the unit
circle this definition (except for the normalization factors associated
with the power spectrum) is precisely the procedure historically used
for estimating the power cepstrum. Thus we may write
Xpc(nT) = (Z(log IX(z) 2)2 (21)
where x (nT) is the power cepstrum.of the sequence x(nT) and X(z)
is the ztransform of x(nT). Consider the signal given below
y(nT) = x(nT) e(nT) (22)
where denotes convolution. The magnitude'squared of the ztransform
of the above equation is then
IY(z) 12 = X(z)12 IE(z)12 (23)
Taking the logarithm of equation (23), we obtain
log lY(z) 2 = log IX(z)12 + log lE(z)(2 (24)
We may now return to the time domain by taking the square of the
inverse ztransform of equation (24) with the contributions due to
the two signals additively combined (provided they occupy different
ranges of quefrency).
pc(nT) = xpc(nT) + e (nT) (25)
Figure 1 shows the computation of the power cepstrum as outlined
above. It is of interest to note that the second squaring
operation will again mix the contributions from the two signals
if they do not occupy distinctly different quefrency ranges (there
will usually be some overlap in practice).
UQJ
0
0J 0
X0 I
4 >
E
4
0,
CL
E
0
"LI
*0
'4
0
41
0
0
O
5
*r
.C
40
0 (
41
0 t.
C)
C 4 
40 r
0
=nO
0E
Ir0
C I
SON
C0 E
=) 0
O O
C r LI..
Z S. U
To further clarify this technique consider a composite signal
consisting of a basic wavelet and a single echo
y(nT) = x(nT) + ax(nTnoT)
(26)
(27)
= x(nT) (6(nT) + a6(nTn T)).
Taking the magnitude squared of the ztransform of equation (27),
we have
jY(z) 2 = X(z)j2 (1+az no
IY(z)!2 = IX(z)I2 I(l+az )I
(28)
Evaluating the above equations on the unit circle (z = ejwT) and
taking the logarithm, we obtain
log lY(e jT)I2 =.log IX(eT) 2 + log(l+a2+2a
= log X(eJiT) 2 + log(l+a2) +
We may now expand the third term on the right
a power series (except for the point values a
to obtain
cos(wn T)) (29)
log(l+ 22a cos(wnoT))
1+a
(210)
hand side of (210) in
= +1 and cos wn T = l)
log (1+ 2a2 cos(wn T)) = (1)m+l (2a cos wn T)m/m
l+a m=1 l+a
(211)
Thus, the logarithm of the magnitude squared of the ztransform of
the composite signal will contain cosinusoidal ripples whose amplitude
and quefrency (that is, the "frequency" of the ripples) are related
to the echo amplitude and delay respectively. If we take the
inverse ztransform of (21TO) and square, we obtain peaks at quefrencies of
10
n T seconds (note, the units of quefrency are seconds, since (210) is a
0
function of frequency the inverse ztransform brings us back into the
time domain) and multiples thereof, These peaks should be detectable
provided loglX(eWT) 2 is approximately quefrency limited to less than
n T, so as not to obscure the peaks. Thus ripples in logjX(ej T) 2
should have a ripple length greater than I/n T. It is now apparent that
the power cepstrum is extremely useful in both echo epoch detection, and
amplitude estimation. It is not obvious from the expansion of (211)
what the relationship is between the echo amplitude a and the heights of
the peaks in the power cepstrum, but it will be seen in the discussion
of the relation between the power and complex cepstra that the
relationship is simple. This will be further discussed in Appendix A.
It should be noted that if we can remove the ripple peaks from the power
cepstrum, it is possible to invert the procedure described to obtain an
estimate of the power spectrum of the basic wavelet, but the basic wavelet
itself can not be recovered since the phase information is discarded
(the power cepstrum is totally independent of the signal phase).
THE COMPLEX CEPSTRUM AND WAVELET RECOVERY
History
The complex cepstrum technique originated as an outgrowth of
the homomorphic system theory set forth by A. Oppenheim.[10] in 1965.
The earlier work of Bogert et al. [8] on the power cepstrum has also
been shown to be a specific application of homomorphic system theory.
The complex cepstrum technique was first presented in R.W. Shafer's
doctoral dissertation [7] in which he derived the technique from
homomorphic system theory., analyzed many aspects of the complex
cepstrum in detail, and used the technique in echo detection, wavelet
recovery and speech analysis. In 1971, Kemerait [3] studied echo
detection and wavelet recovery utilizing the complex cepstrum in
the presence of additive noise, and in the case of an echo distorted
by truncation.
Introduction
The primary difference between the power and complex cepstrum
techniques is that the complex cepstrum technique retains the phase
information of the composite signal. Thus the complex cepstrum
technique may be utilized not only for echo detection, but also in
wavelet recovery.
Formally we define the complex cepstrum as the inverse ztransform
of the (complex) logarithm of the ztransform of a given input
sequence
n
X(z) = Hx(nT)zn (212)
n=
X(z) = log X(z) (213)
x(nT) 2 c log(X(z))zn" dz (214)
where c lies within an annular region in which log X(z) has been
defined as single valued and analytic.
Let us now consider application of the complex cepstrum to a
composite signal formed by convolving two simpler signals
y(nT) = x(nT) e(nT) .
(215)
Ztransforming (215), we obtain the familiar result
Y(z) = X(z)E(z) (216)
Now taking the logarithm of both sides of (216), we have
Y(z) = log Y(z) = log X(z) + log E(z). (217)
Observe that the contributions from the two signals are now
additively combined (just as for the power cepstrum), and that the
phase information of the original signal is retained. We can now
use familiar time domain techniques to achieve separation of the
signals [3]. To complete the calculation of the complex cepstrum
we inverse ztransform equation (217) and obtain
9(nT) = i(nT) + e(nT) (218)
It should be noted that computation of the complex cepstrum is an
invertible operation; thus if the log E(z) contribution can be
"filtered" from Y(z), so that the complex cepstrum becomes
yR(nT) = x(nT) (219)
then we may ztransform (219), exponentiate the result and inverse
ztransform to obtain the signal x(nT). Figure 1 shows the overall
system for wavelet recovery. Note that the encircled operation may
be replaced by a digital filtering operation on the log spectrum treating
it as if it were a time domain function. In fact if this "filtering"
operation can be accomplished by multiplying the complex cepstrum
by some function of (nT) (and ztransforming the result), we may call
it frequency invariant (by analogy with timeinvariant systems),
since the output of the filter is just the convolution of its input
with the "impulse" response.
It is of interest to note that filtering the original signal
y(nT) is accomplished in the quefrency domain by adding the complex
cepstrum of the impulse response of the filter to the complex cepstrum
of y(nT). Also note that filters which are ordinarily unrealizable
in the time domain that is, the impulse response is nonzero for
negative time) are implementable in this manner.
One problem present in the computation of the complex cepstrum
arises from the fact that the complex logarithm is a multivalued
function. If we compute the imaginary part of the logarithm modulo
2n, that is, evaluated at its principal value then discontinuities
appear in the phase curve. This is clearly not allowable since the
log (X(z)) is the ztransform of x(nT) and thus must be .analytic in
some annular region of the zplane. This problem may be rectified
by making the following observations:
(1) The imaginary part of log X(z) must be continuous, and
periodic (evaluated on the unit circle) as a function of w with period
2 since it is the ztransform of x(nT).
(2) Since it is required that the complex cepstrum of a real
function be real, it is required that the imaginary part of
log (X(z)) be an odd function of w.
Subject to the above conditions (and provided the phase curve
is sampled at a sufficient rate) we may compute an unwrapped phase
curve having the correct properties with the following algorithm.
Consider the phase curve (modulo 2r) shown in Figure 2(a). If we
* *
* *
*. a
Figure 2 Phase Unwrapping. (a) phase modulo 2w,
(b) c(k), the correction sequence, and
(c) unwrapped phase curve.
_ _
* 0
are sampling the phase at a rate sufficient to assure that it never
changes by more than a between samples then the.phase may be unwrapped
by adding the correction sequence C(k) to the phase modulo 2T
sequence P(k) where C(k) is determined as follows
C(0) = 0
C(k) = C(kl) 2T If P(k)P(kl) > T
C(k) = C(kl) + 27T If P(kl)P(k) > T
Figure 2(b) shows the correction sequence, and Figure 2(c) shows the
unwrapped phase curve. Alternatively the phase may be unwrapped by
computing the relative phase between adjacent samples of the spectrum.
These phases may then be added to achieve a cumulative (unwrapped)
phase for each point. Both methods have the drawback that the
computation must be done sequentially, i.e., the phase at each point
must be computed before the phase at the next point can be computed.
It should also be noted that if the phase never changes by more
than TT/2 between samples, the phase module T could be computed and
unwrapped with algorithms similar to the above. This is interesting
since it is slightly easier to calculate the phase modulo i than the
phase modulo 27 (the arctangent algorithm is simpler), and many
signals of interest have this property (though noise does not).
There is one class of signals for which phase unwrapping is
unnecessary, namely minimum phase signals. A minimum phase sequence
is a sequence whose ztransform has no poles or zeroes outside the
unit circle. Thus from equation (214)
x(nT) = 0 n < 0 .
(220)
Let e(nT) denote the even part of R(nT). Then
xe(nT) = (x(nT) + x(nT)) (221)
x(nT) = 2u(nT) Re(nT) (222)
where u(nT) = 1 n > 0
= n =
=0 n < 0
but
x (nT) = Z (log IX(z)) .
and thus x(nT) can be computed from a transform involving only the
real part of log X(z). The complex cepstrum determined in this
manner is equivalent to reconstructing the phase associated with the
magnitude of the spectrum using the Hilbert transform. We also see
that for n 0 the complex cepstrum is identical to the power cepstrum
in this case (except for a factor of 2 and the squaring operation).
From equation (220), we see that the complex cepstrum of a minimum
phase sequence is zero at negative quefrencies. Analogously a
maximum phase sequence may be defined (the ztransform has no poles
or zeroes inside the unit circle) which is zero at positive quefrencies.
As will be seen shortly the impulse trains generated in the complex
cepstrum by the presence of a single additive echo are nonzero only
on one side of the origin, and thus will often be referred to as
minimum (or maximum) phase impulse trains. To conclude this section
an example of the application of the complex cepstrum to the single
additive echo case is presented.
In equation (215), let e(nT) = 6(nT)+a6(nTnoT).
That is
y(nT) = x(nT) e(nT) = x(nT)+ax(nTn T) (223)
Taking the ztransform and evaluating it on the unit circle, we have
Y(ejT) = X(ejWT) (l+aejno) .T (224)
Taking the logarithm of both sides of (224), we obtain
?(ejT) = log (Y(ejWT))= log (X(eJ T)) + log (l+aejn oT). (225)
If a < 1 we may expand the right most term in (225) in a power series;
thus
2 3
Y(ejWT) = log X(e WT) + aejnoT e2jno + 3 e3jn .
(226)
Inverse ztransforming (226) we find the complex cepstrum
2
y(nT) = x(nT) + a6(nTn T)  6(nT2nT)...... (227)
Thus the complex cepstrum of the composite signal consists of the
complex cepstrum of the basic wavelet plus a train of 6 functions
located at positive quefrencies at the echo delay (and its multiples)
whose heights are simply related to the echo amplitude. Shafer [7]
has shown that these 6 functions can be effectively removed by
multiplying the.complex cepstrum by a sequence.which is unity.every
where except at multiples of n T, and zero at these points. After.
this comb filter has been applied the complex cepstrum is smoothed
at the zeroed points by averaging adjacent points. The filtered
complex cepstrum is then used to recover the basic wavelet by inverting
the operations used to compute the.complex cepstrum. Had the echo
amplitude (a) been greater than 1, (225) could have been rewritten
y(ejWT) = log (aejnomTX(ejT)) + log (1+ 1 eJnoT) (228)
a
which may be expanded
nT T + T eJnolT j2wn0T
V(ejT) = log (aejno X(ejT 1 T) ej2 (229)
2a
Thus the complex cepstrum will again have peaks at the echo delay
(and multiples),but these peaks occur at negative rather than
positive quefrencies and their amplitudes will be related to 1 rather
a
than a. If these peaks are filtered out and the wavelet recovery
procedure is followed we note that the echo rather than the basic
wavelet is recovered.
From the above discussions, we note that the peaks of the impulse
train in the complex cepstrum may never have an amplitude of greater
than unity regardless of the value of (a). It is also interesting
to note that multiplying the original signal by a scale factor only
changes the n=O term of the complex cepstrum, since the scale factor
appears as a shift in the mean of the.log spectrum. Thus the complex
cepstrum is virtually independent of the signal power.
The Relation Between the Complex and Power Cepstrum
It is obvious that the complex and power cepstrum are closely
related; however, the exact nature of this relationship has never
been fully exploited. This is primarily due to the differing
definitions for the two techniques. From equation (21) the
relationship between the two cepstra becomes apparent,
x p(nT) = (Z(log X(z) X*(z))2 (230)
= (Z(log X(z) + log X*(z)))2 (231)
Assuming x(nT) is real and evaluating its. ztransform on the unit
circle, we find X (z) = X(z ), thus we may write
Xpc(nT) =(2 log X(z) dz 2
X (nT) ( Jlog X(z) zn 1dz + 1 flog X(z1) znl dz)2 (232)
Letting z' = z, we obtain
x (nT) = ( flog X(z) znldz + 2Lflog X(z') ndz2 (233)
pc 21j zj' ) (233)
but the complex cepstrum is by definition
R(nT) =27j log X(z) zndz
Thus we may write (233) as
x p(nT) = ((nT))+ x(nT)2 .(234)
The power cepstrum is just 4 times the square of the even part of
the complex cepstrum. This also follows from the fact that the power
cepstrum is the square of the inverse transform of twice .the real
part of the log spectrum. Equation (234) is of interest since as
pointed out in [3] the power cepstrum is often superior to the
complex cepstrum for epoch detection. Thus in a wavelet recovery
system we may wish to compute both the power and complex cepstrum.
Rather than using the straightforward method outlined in the section
_ __ __ __
on the power cepstrum, we can compute the complex cepstrum, and then
use the simple relation given by (234). This allows the use of one
less FFT in the computation process.
The Phase Cepstrum
The inverse transform of the log phase yields peaks at multiples of
the echo epoch in much the same way that the inverse transform of
the log magnitude does. To see this, let us once again examine the
log spectrum for the single additive echo case given in equation (225)
A T T ~jwn T
Y(eJT) = log (X(e3JT)) + log (l+ae ) (225)
= logJX(eJT) + j Phase[X(ejwT)] (235)
asin n T
+ 2 log (l+a +2a cos wnoT) + jatan +acos0nT .
0
The 4th term in equation (235) produces ripples in the log phase,
just as the third term produces ripples in the log magnitude. Since
Y(ejWT) is the transform of a real sequence, its magnitude (ReY(ejwT))
is an even function of w, and its phase (ImY(ej T)) is an odd function
of w. Thus the inverse transform of Re(Y(ejaT)) will yield the even
portion of the complex cepstrum and the inverse transform of
jIm(Y(e Ja)) will produce the odd portion of the complex cepstrum.
Since the inverse transform of the term log (1+ae j ) produces
peaks on one side of the origin only, the peaks produced by its
real and imaginary parts must be equal in magnitude and must cancel
on one side of the origin while reinforcing on the other (according
to whether (a) is greater or less than one).
II
From the above argument it is obvious that a quantity called the
phase cepstrum may be defined which should be of use in the estimation
of echo epoch times, and amplitudes. Formally the phase cepstrum
may be defined as the magnitude squared of the inverse ztransform of
twice the phase of the ztransform (or equivalently the imaginary part
of the logarithm of the ztransform) of a given input sequence, which
may be written
Xphc(nT) = jZ (2 log X(z)2 log IX(z)1)12' (236)
where the factor of 2 has been introduced to eliminate any normalization
factors in the relation between the phase and complex cepstrum. It
is easy to show from (236) that the phase cepstrum is related to
the complex cepstrum as below
x p(nT) = (x(nT)x(nT))2 (237)
Thus the phase cepstrum is to the log phase just what the power
cepstrum is to the log magnitude. Empirically it has been determined
that the phase cepstrum is less useful than the power cepstrum for
the determination of echo epoch times because of its greater sen
sitivity to noise, but it has proven to be a valuable guide in the
determination of the effects of noise on the signal phase, and may
find application in problems of signal identification. Computation
of the phase cepstrum alone is almost as difficult as that of the
complex cepstrum, since both require the phase unwrapping procedure.
SUMMARY
A new cepstrum technique, the phase cepstrum, has been introduced
which should be useful in the detection of echo epoch times. The
power and phase cepstra have been shown to be (except for a constant)
the square of the even and odd parts of the complex cepstrum.
Because of this close relationship between these techniques, only
the effects of the data handling procedures on the complex cepstrum
are discussed in Chapter IV. The extension of the results to the
power and phase cepstra is obvious. Since the effect of noise on
the phase cepstrum has not been examined previously, this technique
is compared with the power and complex cepstra (in the presence of
additive noise) for the detection of echo epoch times and the esti
mation of echo amplitudes in Appendix A. Some insight is thereby
gained into the results of previous authors.
___
CHAPTER III
THE PITFALLS OF CEPSTRUM COMPUTATION
Since the computation of the complex cepstrum and subsequent
wavelet recovery involve a number of discrete Fourier transforms
together with nonlinear operations, it is not surprising to find
problems associated with these computations not present in ordinary
spectral analysis. This chapter is devoted to the study of these
problems,some of which have been noted by other authors; others
were previously unknown.
LINEAR PHASE TERMS
The phase of the ztransform of a composite signal will often
contain a linear trend. This may be a property of the signal under
analysis, or may be due to the signal being delayed (that is, not
starting at the beginning of the record). In such cases the ztransform
of a typical input signal can be rewritten
X(z) = zrX'(z) (31)
where zr is the term which contributes the linear part of the phase
curve. Then
X(z) = log X(z) = log zr + log X'(z) (32)
Consider the contribution of the term log z to the complex cepstrum.
XLp(nT) = 2j flog zrznd (33)
Evaluating the above expression on the unit circle z = e T, we
obtain
Lp(nT) = (jwTr)ejnWTjd(wT) (34)
Integrating, we find that the contribution of the linear phase term
to the complex cepstrum is
xLP(nT) = 0 n=0 (35)
S cos Tn ntO
n
The linear phase component thus causes a rapid oscillation in the
complex cepstrum on which the wavelet and echo information rides.
If the linear phase component is sufficiently large, it will dominate
the complex (and phase) cepstrum,masking important characteristics.
Figure 3 shows the mean square error (MSE) of the recovered
wavelet as a function of the delay in a composite signal. This delay
only affects the phase of the composite signal (introducing a linear
component). Note the rapid increase in MSE as a function of the
delay. For comparison the MSE of the recovered wavelet is shown
when the linear phase term is removed prior to computing the complex
cepstrum. No explanation is available for the very low MSE observed
at the odd delay times. For the composite signal considered (echo
amplitude = .5) the echo peak becomes undetectable in both the phase
and complex cepstra at a delay of 30 sample times (when no linear
* Complete Phase Removal
* No Phase Removal
SI I I I I
15 30 45 60 75
Signal Delay
Figure 3 Effect of a Delay in the Composite Signal on
the MSE of the Recovered Wavelet.
phase removal is performed). Clearly the removal of linear phase
terms is required for an effective wavelet recovery system.
Linear phase removal results in a dramatic change in the
appearance of the complex cepstrum and also displaces the recovered
wavelet in time (unless corrected for). The linear phase component
is removed by calculating the unwrapped phase of the N/21 point,
using this to compute the slope of the linear phase term, and then
subtracting the linear phase term from the signal phase. For
the purpose of linear phase removal the unwrapped phase of the
N/21 point can be truncated to an integer multiple of r [3].
This allows the correction of the recovered wavelet (for linear
phase removal) by shifting it an integer number of sample times.
This generally still leaves a small linear phase term present in
the log spectrum. Figure 7 shows the complex cepstrum for this
type of linear phase removal. Note the effects of the remaining
linear phase term are clearly visible.
Figure 14 shows the complex cepstrum of the same signal as
Figure 7, but here the linear phase term was removed entirely prior
to computation of the complex cepstrum. Note the radically different
appearance of the complex cepstrum. In this case, the linear phase
term must be reinserted in the log spectrum of the recovered wavelet
during the inverse complex cepstrum operations.
It should be noted that if the linear phase term is not entirely
removed from the log spectrum of the composite signal then the
echo peaks cannot be smoothed (without introducing considerable
error into the recovered wavelet) by averaging the points adjacent
to them since these points have contributions from the linear phase
Example
Incomplete Linear Phase Removal
Figures 4 through 10
This group of figures illustrates the computation of the cepstra,
and wavelet recovery when the linear phase term removed from the
unwrapped phase curve is based on the unwrapped phase of the N/21
point truncated to an integer multiple of r.
The composite signal (256 points) is
y(nT) = x(nT) + .5x(nT30T)
where x(nT) = nT enT 0 5 n < 64
Figure Number Figure Title
4 Composite Signal
SNR = 40 dB
5 Unwrapped Phase Curve
6 Log Magnitude
7 Complex Cepstrum
8 Phase Cepstrum
9 Power Cepstrum
10 Recovered Wavelet
MSE = 7.97 x 107
28
I   
I I
M
I
rrrll~l~lrr~r~ ;F
o
I 3
Si i a
: 
*I 4. i
I I I I I *
C)
. 
I ao i
a l   /

3 Ca o 34 a'
,rrrr:~clr*: r~
I * *
* a I
* S. S S
S. 5 5 r S I*1I II* J.u*1
* S
* 5 a S i
S S
* SI S
*I a
S
"
ft *
S a
I
: *( xt r n n ^LII i  ~r *.1in _ ^__ r r  __. ^___..,___
* ma : :" i
* S
i
* S a S
ft
a S S
S
fC__ II"I.rLeDC *_
*
*~ ~ a, S
S
tf
__ .. : .
s. S a a K g
S. *
S S
5 S. 5
*' a : S
* a S S
* a
*l S
* S
*> ...... .. . __________
~ ~ ~ a *S
*~~i *5 5
* S C 5 5 5 5; ~ aU~
*r S
*
~~r~rrr~
C
,,~.
I
I
a a
a == z a: i.
a a
S a
* a
* a
* a
I a
* a
I S
4 a
* a
* ..a ~. .. .,
*r a a
* a a a
a a p 3
a a a a a
*~r a .
* a a a. a a
~ ~ l __ CI 
o) I a
a a
* .
* ~ ~ a *. a
* ? *c .e.
* a a a
* a g a a
* a a
C a.a
* a
ac a
* a. a
*~~ ~ a
*~cZ L I a ~  a a a a
* a
* aaa
*i a. a
a a
a a
" 01
I ?
0
0 0
o
0 .
.?
L
0 0
o
U C    ^ ^^. s
* a a a
a
a
S~~ I a
* ** a a.^^ m^ c ill~B *, 
* a S
*^ M a a ^ a a 
,.* a b* S
* a S
*. aSa
I a a aa
a a.
a a a S.
a
a0 a
I
~ i
i '
*5a a
fta as a
a a a
* *^_ a s
Sea a a
S. a p
a. ~ a
a a a a
a at a*rrrr~rrrrorr
a a a
a. a 4
a a
4 . .re...e
* a j a m
4,
S S
S S
S S
')IP~.........) ...........
: i;
.. I .
* 9
i 5
S O5
0 3 5
0 0 
* 0 0**0*55O**0*C S ~00~c0*0r0*~~ 4 r4 ..r...e
* S 5 5 00
S 5 3 0 5
S S 5 0 0
*"LLCrC~rrC))IC) 3 5 8 0 5
~ ~ 3 5 0
S 5 0
~ ~ S S 5 0
5 3 5 0 0
*rr~u crr ~c~r~r S 5 5 0 0
* . .... o o.aS. oe. I *4 f .5 4 0
SS 5 0 0
S S 5 0
5 3 5 0 5
S 5 5 0 0
*'~~t~~~ Sar~rr 5 5 0 58 r
* S S 0 O
* S 0 0
S 0 0
~ ~ ~ 3 0
*I S 5
S 5 5 SO
S 5 0 0
* SICIC r~~+,, S V 0 0
~ ~ ~ 5 5 0
* S S
* S0
5 5 S SO
~ ~ 5 S 
+5 0 n 4 5 0 a.. S r...I. 0 05I ~ .~
* S 5 5 S 05P~r CCF~U; ~C
* SS 0
05
* S 5S 0
S S 5 00
5 3 5 0
5 .4.0?r I 00... S
0 0
.. o *a..0 00 50
* 50
Oo"
0s
SC
S.
5.
*I.
5 3 0 *40 *
S SS 0 S
~ ~ ~ S S
~ ~ ~ S *
*~~ S S o S
~ ~ ~ S S 5
S 0 55 0
*r 5 0 5 05
S S 0 5
o a a a
0
SE
S
x
u
I'
3
SL
n
o
. .. S S .. S .
I
j
U.
:i
__ ^ __ 1 "
U U
I S: U U
I U I S I4
:a
IU
.
I S I U U ft
i S U
: I U
ii
I U U S 4
I U U a
U
I U
I S U S
I S
U II U
*ur c~~r S~ C~ u~ l S~ IIFW I U CL~C II U 4~CI
*~ I t
U U I U
SI U
I I U 4
I I 4
U I S 3 4
I U U I S
I.I I I 5: 4
*^ U S U U
I U I
hS
II 4
* I. I S
*IUII 44
I SII4S4
I I S4
I S US
I I U S4
I S II US
* I I 4
*. *< I* 5' 5
! t
i i j i _i
i '
* ,
i
i i
* *
e ii" l "" ~
I i
, ^rCI ^*nr^r,' ^.^. ^I ^ 'Y IO.
1 t1
: 
_*j _ a_ a a a
*
a aa
a
a a a a a a ^ .
Q C
D C
O O
0
O
O O
O 5
4f ~ fl S
S S S
S SS S
S S S S
~ ~ S
S S S
* S S
* S S5 5
S S S
S S S S
*~ S
S SS S
* S S
S SS S
* S S
S S S S 5 0Lu~~)~)+~L
S~ S S S
* S S S S
S SS S
S SS S
* 55 S I 4.O.B4ruu ~ 0 B.B45i
* S S S
S SS S
S B5 5
*r~,~r c c~~c~rr S S S S Bu~
* S S
* S: S S S
* S S S
*~uo S c~~rrr~lr S S r~rur~~ S
* S S S S
*~ S S S
*
*~~ S S
*~ S S
*~ S S
S S S
* ~;. S S S
* S
* S S
* S SS S
~ ~ S: S
*rsoLp~ro~ICI~ S
* S S S S
* S S S
* S B S
S 5 5B,
S S S S
*~~ S S
* S S SS
* S S SS
~~~ S SS
* S S S S
~~rr ~ cr  ,...n s., S S
* S S SY S
*N S SS
* S
S
S SS S
41
1
a r
>I
o X
0
* x
00
S
:3
tm LI.J
*r
Example
Complete Linear Phase Removal
Figures 11 through 17
This group of figures illustrates the computation of the cepstra,
and wavelet recovery when the linear phase term removed from the
unwrapped phase curve is based on the unwrapped phase of the N/21 point
(no truncation).
The composite signal (256 points) is
y(nT) = x(nT) + .5x(nT30T)
where x(nT) = nT enT 0 5 n < 64
Figure Number Figure Title
11 Composite Signal
SNR = 40 dB
12 Unwrapped Phase Curve
13 Log Magnitude
14 Complex Cepstrum
15 Phase Cepstrum
16 Power Cepstrum
17 Recovered Wavelet
MSE = 8.01 x 107
0
I I I I j
Itt
I i
I I
J i !I
I I IP
a ac
i i "
* 
I I s
I I
4.tM b * k44JJl4 J* 4* UHJ kMlr
Sa i i
Sft^  ^ t  .. ^ .... ..  ft .
I I I a
a a a
a I I I 
a an
a 0)
I i I *
aIrlr,,rr,+rrrr~~r S
a~~ a aa
a a a a *. o aIo
a a a
I a aa a .a Vi ILL/
a a aa a tft
I t ft 4* f
a a a a a
a~~~ a t
a a a
C C 0
0 0 0ad 0
Pa (4 d 0
0r 0 0 0 0 0
I ;
I I I. I I
IIi I I
I I I j a
. I
i I I I IC
I I I I
a a a as
I i i i
I I I I
i I I I I
I I a
I I0.
SI I I
cii
S* I I 1
" i i io
.   
SI i
Ci
c
I I

a C I I I.
*s,
1 1 I I *I
(* I I I I1
a a. I I I
I I C1I I C
I I a i a.
r  c  c  
m o o
0
6 j ( 
I 1 0
SI ":I I *
! i r i 
I I o
~ ~ t : .
' ^ .^ i......... .... I. .. ^ ^
I 1 I o I
I I I i 6
I 6 6 6 i6 s
^._.^      ^0.0
SI !.
i J I 6 * o
*I * 60
*i
I I i I o
I I I I
~
Ioi
i i i E
0 6
I i e 1 i
i *' I lo
I I *I
I... 1 1
i 6 I i 6 6 i .
1 i I I~ I I o
6 6I 6* 6
6 I 11 0
1* 0 I!
I 6 I I l
I 6
I I  1 1
i~~ *i i ri
6 6 f 0. 0 60
*0.
: ~ I
I 60
6 0,
It **i I
CI*
O P1
6 6
39
. btMJUf ^ >MlM*h4M k  .M  h 
I :" Io
4..
SI I H
I iI I I
O *
^  
  . . .  .. ....... .
4M*1l1.11M 1,  1 f    
I C
I I I IC
j
I E
1 I
I. 1o ,
io D
I I 0
* i U
*n o o  4" o l    
  I o 
I I iJ S
i j i I 4 .
0 ' C N N
5
j!
* 4
jaa
6 4,
.S
0
0
a a
40
0
ftj
0
I a
0
a0
S a
a a
aC aU(
a c
I a
a *
I a 61
a"rQ CCI aC~Y I a IC((CCCL~ ~ l aC a
a a a a a
a a aa a
* a a a a
a a a a
* 5 .65
* *" C~+r~I *fl~r C ~ ICa. a..*.a).,*.
a a a S .6
*~ ~~ 6
*an
I ~ ~ a
a a a
a
*I a a
I~ 4
aCC~~"~+rCr~ 6 a a a aIC ~ L~C~)(~I
0^
i i
Sa
I I
I~ a
i
a a
a a aU, Ur ~ CCI~IUCUC~lC
a a a
a a a
a a a
a a
I
a~L~'l' aCCC I a aCCtUj C"~CCI
CI C C C
P
R
C,
F
O
N
O
O
Fl
O
O
P
O
D
N
Y1
I
~
41
o
I s
* S__ a_ S S a 
_____ ___ __ __ .,
I
\ i
r ___CI~ U I CI~ ____ __ __ __^ <
* a0
; :
I 0.
SS I a a
* 0\ t .0 0 S 0000 S
sI
S
a a*
SI '
I ;: __ ^
_~ 1 ._ ____s
* aa
* S a
I : M ^
a
a a a
I S S S*
* a a a a a 0
* a a a
* a a a a 5
* s
* a SS
~~~ a
S S a a
* S :
*~ a a aa
* a a a .a a a
*,,r ar~rre au ,w~~r~~ a a a at'
a a a
* a a a a 
SL S S a
*~ a a a
  i .
SU :
I 
I I t
4.1llltl*ill '   
 t11 U kr*f~fjJ Ml44jJh~' f~cMr~* <<. * *
I I 1 1 *
 a  *
i"
i. o
ii o
I i I ,I
I 
I   a. ..
a a' a <* f
i,,,,
a a a a a
component which are opposite in sign to the contribution of the
point being smoothed. To smooth the nth point properly when a
linear phase term is present the contributions of the n+2 and n2
points must be averaged. This form of smoothing results in a
somewhat smaller MSE (than averaging the points adjacent to the
echo peaks) even when the linear phase term is completely removed.
Subject to the above smoothing considerations it has been
determined that the MSE of the recovered wavelet is not significantly
different for the two forms of linear phase removal, though as
mentioned previously the appearance of the complex (and phase)
cepstrum is quite different. This undoubtedly will affect the
detection of small amplitude echo peaks and echo peaks occurring
at low quefrencies. Thus the second form of linear phase removal
is preferred, and is employed in the computer experiments discussed
in this dissertation.
PHASE UNWRAPPING ERRORS
Previous authors have tacitly assumed that the phase of the DFT
of a given sequence may be unwrapped unambiguously. Only recently
has it been pointed out that errors may occur in the phase unwrapping
process [5]. Since the arctangent routine used in the phase
unwrapping algorithm computes the phase modulo 2n, it is evident
that phase unwrapping (explained in Chapter II) is performed correctly
only if the change in phase between samples of the ztransform of the
given sequence is less than Tr. The DFT is just the ztransform
(evaluated on the unit circle) of the sequence sampled at values
w = na =2r/TR where TR = NT is the record length under consideration.
The sampling theorem (in the frequency domain) indicates that Aw
must be less than or equal to 27r/TR in order to reconstruct X(z)
exactly from its samples. Thus the DFT samples the ztransform
at the minimum rate necessary for reconstruction. Sampling at this
rate does not, however, assure us that the phase of the ztransform
(evaluated on the unit circle z=ejWT) does not change by more than
a between samples as can be seen by the following simple examples.
Example 1
Let x(nT) = 6((nN/2)T) + 6((nN)T) (36)
NT
X(z)= X(e T) = Je + ejNT (37)
z=ejWT
Consider the phase diagram of X(ej T) shown in Figure 18. It is
evident from Figure 18 that even if the ztransform is sampled at
twice the minimum rate (A = < ) the change in phase can be
R R
greater than w. Thus the phase unwrapping algorithm will yield
incorrect results. To see the effects of this on the computed
phase curve consider Figure 19. Obviously the computed phase curve
differs considerably from the true phase curve.
Example 2
Let x(nT) = y(nTnoT) [O,N) (38)
= 0 otherwise
jwn T c.T
X(ejwT)= e o yj(eJ ) (39)
As expected the phase of X(e j) is the sum of a linear phase component
and the phase of Y(ejT ). If w is sampled at the minimum rate, that
n2s
is, w = nAwn = n then
(b)
Figure 18 Phase of X(eJWT).
(a) w=O, and (b) w=T/TR.
671
Figure 19 Phase Unwrapping Errors.
(b) phase modulo 2i, and
I
(a) true phase,
(c) unwrapped phase.
jn2r
Xs(eJ) = X(e ) = ej2n(no/N)y (eJ) (310)
Thus if no > N/2 the linear component of the phase will change by
more than 7 between samples and unless the phase of Y(ejWT)
counteracts this change, the phase unwrapping algorithm will yield
erroneous results. This was observed in computer experiments when
the composite signal is delayed by more than half the record
length. As expected this not only reduces the echo detectability
in both the phase and complex cepstra, but also severely distorts
the recovered wavelet.
Fortunately,provided the composite signal starts relatively
early in the record, many signals of interest can be unwrapped
correctly since their phase changes slowly. Furthermore many signals
only occupy a portion of the total record, and thus the DFT yields
samples of the ztransform spaced more closely than dictated by the
minimum sampling considerations. Even in the absence of the above
conditions, this problem can be alleviated by extending the data
record with zeroes as this increases the sampling rate of the
ztransform (this will be shown in the next chapter).
These errors often occur when unwrapping the phase of a composite
signal in additive noise at low signaltonoise ratios (SNR), since
the noise may dominate portions of the spectrum and the phase of
the noise spectrum will change rapidly from sample to sample. This
error in phase unwrapping seems to be actually beneficial in this
case in that it always limits the jumps in phase (to less than 1)
when in actuality the jumps may be somewhat larger. This topic
will be discussed in more detail in the section on the addition of
zeroes in the following chapter.
ALIASING
Even though the DFT yields adequately spaced samples of the
ztransform (X(z)), it is not surprising that aliasing may be a
problem in computation of the cepstra, since the nonlinear complex
logarithm operation will introduce harmonics into X(z).
Consider, for example, the complex cepstrum. Because of the
harmonics introduced into X(z), the complex cepstrum will be infinite
in extent. The DFT will yield a finite N point sequence, thus the
NT NT
contributions from quefrencies nT > and nT < will alias
into the interval ( N n).
The aliasing of the cepstra of the basic wavelets considered
did not prove to be a problem. This might have been expected since
their time duration was never more than 25% of the total record
length.
One manifestation of aliasing is an ambiguity in the determination
of the echo epoch time. A peak at n T in the complex cepstrum may
be caused either by an echo delayed by (Nno)T or by an echo delayed
by n T. This ambiguity cannot be resolved unless the relative echo
amplitude is known to be either greater or less than unity (unless
the original data record is extended by the addition of zeroes).
Aliasing of the echo impulse train can also be a problem. The
aliasing of the higher order terms in the log expansion given in
equation (225) is often observed. This is especially troublesome
if (a) is close to one since the amplitudes of the peaks in the
impulse train do not fall rapidly in this case.
Aliasing is also evident when additive noise is introduced into
the observed record. Since noise is present throughout the observed
record it is not surprising that it produces contributions at high
quefrencies in the complex cepstrum and thus may cause aliasing.
In the noise analysis of Appendix B, it is shown that aliasing
is the primary mechanism through which noise is introduced into
the complex cepstrum.at negative quefrencies for high SNR.
As will be seen in the next chapter,adding zeroes to the input
data sequence reduces the aliasing of the complex cepstrum. In
fact, if the number of zeroes added is greater than or equal to
the original record length the ambiguity discussed above can never
occur since the echo delay will always be less than one half the
augmented record length.
Aliasing can also be reduced (for a given composite signal) by
choosing the record length NT as large as possible subject to the
constraints on the number of points analyzable and on the minimum
sampling rate. That this is the case becomes evident if we recall
I 1
that the.samples of the log spectrum are spaced Af = T =T apart;
thus increasing NT increases the "sampling rate" of the log spectrum,
and should minimize aliasing in its inverse DFT. This gain is offset
somewhat by the fact that increasing the duration in time of the
observed data increases the "long time" (high quefrency) contributions
to the cepstra from noise (assuming the composite signal remains
constant). Furthermore the shorter the interval of the data record
that the composite signal occupies,the lower the SNR for a given
noise environment.
OVERSAMPLING
One problem caused by additive noise occurs when the composite
signal is oversampled. Outside the signal band (although.these
components may have been greatly attenuated prior to sampling) noise
often dominates the spectrum. This presents no problem in ordinary
spectral analysis as these components contain little power, but
they may have a marked effect on the cepstra (since the complex
logarithm of the spectrum is taken). Because of this nonlinearity
the regions of low power may contribute as much or more to the
.cepstra as the signal band, and therefore will degrade both echo
detectability and eventual wavelet recovery. This effect is especially
evident in the phase cepstrum since the phase of the spectrum outside
the signal band may change rapidly from point to point. It is difficult
to assess the effect of oversampling on the recovered wavelet, but.
the appearances of the complex, phase and power cepstra are observably
degraded by oversampling.
Oversampling also aggravates the phase unwrapping and aliasing
problems previously noted, since it shortens the data record (if the
number of points is fixed) which in turn implies that the samples of
the log spectrum are spaced farther apart.
SUMMARY
Four problems connected with cepstral analysis have been examined.
These problems are linear phase terms, phase unwrapping errors, aliasing
of the cepstra, and oversampling.
__
50
The presence of a linear phase.term distorts the appearance of
the complex and phase cepstra even when partially removed as in
reference [3]. Such a term degrades not only echo epoch detectability
in the phase and complex cepstra, but also wavelet recovery. The
complete removal of linear phase contributions to the cepstra is a
necessity for satisfactory wavelet recovery.
Phase unwrapping errors are observed to occur if the log phase
changes by more than. between samples. This can be quite detri
mental to echo epoch detectability in the complex and phase cepstra,
and to wavelet recovery.
Some aliasing of the cepstra is always present since the cepstra
are derived from samples of the log spectrum. This results in an
ambiguity in the echo epoch determination. Both the aliasing and
phase unwrapping problems can be alleviated by extending the original
data record with zeroes.
Finally, oversampling the input data record is observed to degrade
wavelet recovery when additive noise is present. Oversampling also
aggravates the phase unwrapping and aliasing problems, and thus
should be avoided.
CHAPTER IV
THE EFFECTS OF SOME DATA PROCESSING TECHNIQUES
ON THE CEPSTRA
Care must be exercised when computing the complex, phase, and
power cepstra. The repeated use of the DFT may raise the problems
of aliasing, leakage, and the picket fence effect. Since the complex,
phase, and power cepstra are nonlinear techniques, it is not clear
that ordinary data processing techniques, normally used to improve
spectral analysis, are applicable. Two such techniques, windowing
(normally used to reduce leakage) and the addition of zeroes (normally
used to reduce the picket fence effect),are examined to determine
the effects of their application at various stages in the echo
detection and wavelet recovery process. Primarily the effects of
these techniques on the complex cepstrum are discussed, since the
extension of the results to the phase and power cepstra is usually
obvious.
WINDOWING THE COMPOSITE SIGNAL
When the input data record is windowed the applicability of
cepstrum techniques is no longer evident. Consider for a single
additive echo
y(nT) = [x(nTk T) + ax(nTk Tn T)]w(nT) (41)
where w(nT) = 0 n > L1 or n < 0
x(nT) = 0 n < 0
w(nT) is the windowing sequence. The basic wavelet x(nT) has the
argument (nk )T, since we wish it to begin at some unspecified
point (k T) because it is not known where (if at all) in the data
record the composite signal begins. Equation (41) is rewritten
y(nT) = [x(nT)*(6(nTkoT)+aS(nTkoTnoT))]w(nT) (42)
Ztransforming (42), we obtain
k n
Y(z) = [z o(X(z)(l+az O))]*W(z) (43)
Thus, the contributions of the basic wavelet x(nT) and the echo delay
cannot be separated by taking the logarithm since the term in brackets
is convolved with W(z). Fortunately many practical situations permit
the extraction of the basic wavelet (though there is some error)
and detection of the echo arrival time.
Experiments have been undertaken to determine the effects of
windowing the composite signal on the wavelet recoverability (as
measured by the MSE between the recovered wavelet and the true
wavelet) and on the echo epoch detectability. These effects are
ascertained in both the noise free and the additive noise case.
Unless otherwise noted the composite signal examined is of the form
y(nT) = x(nTk T) + ax(nTk TnoT) 0 n [ L1 (44)
and x(nT) =nTebnT 0 n < 64 (45)
x(nT)= nTe"bnT 0 n < 64 (45)
with L = 256
a =.5
n = 30
k= 0
b = .1, .5, 1.0
T = .333 .
A typical experiment is conducted as follows: a composite signal as
given above is generated. Bandlimited random noise is added to this
sequence and it is windowed. The complex, phase, and power cepstra
are computed. Peaks due to the echo train are smoothed in the complex
cepstrum, which is then inverse transformed to obtain an estimate of
the basic wavelet. This estimate is corrected for the windowing by
multiplication by (w(nT))l where w(nT) is the windowing sequence.
The mean square error (MSE) is computed between the recovered and the
initial basic wavelet over the entire record length.
The Exponential Window
The exponential window
w(nT) = an 0 S n < L
a < I
(46)
= 0 otherwise
has been proposed by Shafer (7) to reduce any error associated with
truncating the echo. To some extent this window also reduces
leakage, but since its properties are quite different from those of
the commonly used windows, it will be studied separately.
Experiments were undertaken to determine the effects of this
window on wavelet recoverability and to determine its effects on
echo detectability both in the noise free and the additive noise
cases.
The following observations are.made on the effects of this
window:
(1) No difficulty is encountered detecting the echo epoch,
though the peaks in the power, phase and complex cepstrum are
reduced in height. The echo detectability threshold in the case of
additive noise is about the same as observed for the rectangular
window (about a signal to noise ratio (SNR) of 8dB).
(2) The MSE of the recovered wavelet is somewhat better than that
obtained with the rectangular window at low SNR (14 dB and below),
but not nearly as good at higher SNR. A comparison of the MSE is
given in Figure 20. At high SNR the MSE in the recovered wavelet
is comparable to that for the rectangularly windowed case if it is
computed only over the known duration of the basic wavelet. Thus
the MSE results shown in Figure 20 reflect the distortions introduced
into the recovered record outside the signal duration by correcting
for the exponential window.
Interpretation of Results
The following derivation is presented in order to clarify the
effects of the exponential window. Consider again the ztransform
of equation (41), it is:
L1
Y(z) = [w(nT)x(nTk T)+aw(nT)x(nTk Tn T)]zn (47)
n=0
n Ln 1
n '0(o)
Y(z) = X (z)+az o x(nTk T)w(nT+n T)zn (48)
n=n
* Rectangular Window
* Exponential Window
40 . 
40
Figure 20 MSE of Recovered Wavelet when the Input Data
Record is Exponentially Windowed
102+
103
20
SNR(dB)
L1
where X (z) = > w(nT)x(nTk T)zn .
n=O
n Ln 1
Y(z) = X (z) + az o x(nTk T)w(nT+n T)zn+E (49)
n=0 o
n 1
where El = az o > x(nTk T)w(nT+n T)z
n=n
0
Now letting w(nT) = an 0 S n < L
= 0 otherwise
n n Lno1
Y(z) = X (z) + az o a o x(nTk T) anzn+E1 (410)
n n L1 noLL no nn
= X (z)+az O[a 0 o x(nTk T)anzzn0 a x(nT+LTn Tk T)aa nz
n=O n=O
+ E (411)
The term denoted El is an error term. A close examination reveals
that the error E arises because it has not been assumed that the
basic wavelet begins within the windowed record. If this assumption
is made (that is, k must be positive and less than L1), then the
error term El vanishes. Henceforth this assumption is made unless
otherwise noted. The second term in brackets also is an error term.
It arises because the echo has been truncated. It is the error due
to this term that this window minimizes (due to the presence of the
factor aL). At present it is assumed that the duration of x(nT) is
less than (Lno o)T so that the truncation error vanishes. With
the above assumption equation (411) becomes
n n n n
Y(z) = Xw(z) + oaz o Xw(z) = X w(z)(1ha Oaz 0) (412)
n n
Y(z) = log Y(z) = log Xw(z) + log (l+a'oaz o) (413)
n
Clearly the presence of the factor a 0 alters the heights of the
peaks which are observed in the complex cepstrum. If a < 1, the
peaks in the complex cepstrum are reduced by the factor (a o)m
where m is the order of the peak (as compared to the rectangularly
windowed case). This is precisely what is observed in the experi
mental section. If a > 1, and a na > 1 then the peaks in the
nom
complex cepstrum are enhanced by the factor (ao ) If a > 1, and
a na < 1, then the peaks may be either reduced or enhanced, and
they occur at positive rather than negative quefrencies. The basic
wavelet rather than the echo is then recovered.
It can be seen from the above discussion that aliasing of the
echo peaks can be reduced by choosing a properly. Another possible
advantage to using this window occurs in the case of multiple echos.
Kemerait [3] has indicated that when the sum of the echo amplitudes
is near one,ambiguous peaks occur in the cepstrum. Again by choosing
a properly this situation can be avoided. Assuming that the contri
n n
bution to the log spectrum from the term log(l+a z ) can be
filtered out (by transforming and smoothing the corresponding peaks
in the complex cepstrum), the inverse wavelet recovery procedure
should, after correcting for the windowing, yield the basic wavelet.
Experimentally this is found to be true over the duration of the basic
wavelet, but the windowing correction appears to introduce some error
into the recovered record outside this interval and thus the MSE
results shown in Figure 20 are somewhat greater than observed for the
rectangular window at highSNR. As stated previously,at low SNR
the MSE is substantially lower than that obtained using the rectangular
window. It is also noted that the echo detectability threshold is
about the same as observed for the rectangular window, even though
the peaks in the complex cepstrum are reduced in amplitude. Both
of the above results are undoubtedly due to the fact that the signal
occurred near the beginning of the record and thus is not reduced
as much by the windowing as is the noise which is spread throughout
the record. Had the signal occurred much later in the record, it
is expected that the MSE would increase considerably, and that the
echo detectability would suffer. This has subsequently been verified
experimentally. Essentially, the'cost of reducing the error due
to truncation with the exponential window is that the portion of the
signal occurring late in the:data record is also discriminated against.
The Common Window
Three windows (Hamming, Hanning, and Tapering) commonly used to
reduce leakage in spectrum analysis are examined to determine the
effects of their application. Since these windows exhibit similar
properties, they are discussed as a group in this section with the
Hamming window being discussed first, and the other windows then
compared with it.
.The following observations are made on the effects of the
Hamming window.
(1) Several signals of the form given in equations (44) and
(45) with (b) ranging from .1 to 1.0 were generated without additive
noise. In all cases the echo epoch is detectable, though the height
of the peaks found in the complex, phase, and power cepstrum are
quite different than are found when the rectangular window.is used.
Substantial peaks at the echo epoch time are often found at both
positive and negative quefrencies in the complex .cepstrum with the
negative quefrency peaks being larger (the rectangularly windowed
sequence has peaks at positive quefrencies only). The peaks are
somewhat smeared. For the echo delay used (n =30), satisfactory
wavelet extraction is not possible except when b=1.0; in this case
extraction of the echo rather than the basic wavelet produces the
least error (since the negative quefrency peaks dominate the complex
cepstrum) though considerable distortion is still evident in the
recovered wavelet.
(2) Since wavelet recovery is only possible when b=1.0, a
number of experiments were conducted with b=1.0 to determine the effects
of the window when noise is added. Figure 21 shows the MSE of the
recovered wavelet as a function of the SNR when the input data
record is Hamming windowed. For reference the MSE obtained with the
rectangular window is also shown. It is noted that the MSE is
sometimes reduced for the Hamming windowed composite signal by
smoothing not only the peaks, but several adjacent points as well.
For this Hamming windowed composite signal the peaks in the
.complex cepstrum occur at negative rather than positive quefrencies.
In the no noise case the negative quefrency peaks are larger than for
the corresponding rectangularly windowed signal. Even with this
increased peak height in the no noise case the SNR at which the echo
epoch becomes undetectable is about 12 dB greater than found for the
rectangularly windowed case.
__
* Rectangular Window
* Hamming Window
10io
 
40
SNR(dB)
Figure 21 MSE of the Recovered Wavelet when
Are Hamming Windowed.
the Input Data
10
Only a few computer computations were made using the Hanning
and Tapering windows, but general trends have been ascertained from
the results. The Hanning window produces considerably worse results
than the Hamming. Likewise, the Tapering window gives poor results
if a major portion of the composite signal lies in the region of
tapering; however, if the composite signal lies in the untapered region
the results are almost identical to those for the rectangular window.
One additional drawback of these two windows is that they are zero
at the endpoints and thus it is impossible to correct the recovered
wavelet at these points. It is noted that correction for the window
when w(nT)<
introduce some distortion.
Interpretation of results
The following theoretical analysis is intended to clarify the
conditions under which wavelet extraction is possible after using
one of the foregoing windows. Consider (47) which is rewritten as
L1 Ll
Y(z) = w(k T)~2 x(nTk T)zn+aw(n T+k T)7 x(nTk Tn T)zn
n=O n=O
L1 Ll
+ 2 x(nTk T)[w(nT)w(k T)z+a x(nTkTn T)[w(nT)w(n T+k T)]zn
n=O n=O 0 0 0 0
(414)
The third and fourth terms of equation (414) can be considered error
terms. It will be seen below that extraction of the wavelet is
possible if these terms are small. This requires the window to be
approximately constant over the duration of x(nT). For the windows
under consideration this requires the duration of x(nT) to be much
less than the record length. The error terms vanish completely when
the window is rectangular. The sum of the third and fourth terms
will be designated as error term E2. Equation (414) is now rewritten
as
n
Y(z) = [w(koT) + az o w(noT+kT)] X(z) + E2 (415)
L1
where X(z) = : x(nTk T)zn
n=O
aw(n T+k T)
If = WoTTo >1, then
n E2
Y(z) = log Y(z) = log[az w(noT+k T)X(z) + E2 ] + log(l+z0o/) .
00 l+z n/
(416)
The term log (l+zno/B) may now be expanded as we have done previously
to obtain
n E2 "o 2no
Y(z) = log[az w(nT+k T)X(z) + l+z + 7 . . (417)
0 1+z /e 2 0
Thus the peaks in the complex cepstrum occur at negative quefrencies.
If these peaks are removed, then the recovered wavelet will be
n E
yRW(nT) = Zl(az w(n T+k T)X(z) + ) (418)
S0 0 l+z /2
1
Expanding n in a series, and noting the origin of the error
l+z /8
term E2, we may rewrite equation (418) as
yRW(nT) = aw(n T+k T)x(nTk TnoT)+[x(nTk T)(w(nT)w(k T))
+ ax(nTk TnoT)(w(nT)W(n T+k T))]*[6(nT) 6(nT+noT)
+ 26(nT+2nrT)...] = w(nT)ax(nTk Tn T)+[x(nTk T)(w(nT)
w(k T))]*[6(nT)B16(nT+n T)+26(nT+2n T)...]
+ [x(nTkTnoT)(w(nT)w(n T+kT))]*[B (nT+n T)
+ 2(nT+2n T)... (419)
After correction for windowing, and a few simple manipulations,
equation (419) becomes
w(koT)
YR(nT) = ax(nTk Tn T)+x(nTkoT)(1 0
1 w(nT+noT) 1 w(noT+koT)
w(nT) + W(nT)
w(nT+noT) w(koT)
x(nT+n Tk T)( w(nT) w( nT
1 w(nT+2noT) 1 w(noT+koT)
w w(nT) w(nT) )+ "' (420)
Thus the recovered wavelet consists of the echo plus a number of
distorted and shifted replicas of the basic wavelet. The above
analysis agrees well with the experimental results. For the case
considered (no=30, a=.5 with the Hamming window), =1.25. Thus peaks
are expected at negative quefrencies and since = .8>a, these peaks
should be larger than for the corresponding rectangularly windowed
case. The echo can be recovered by filtering these peaks from the
complex cepstrum. The distortions suggested by equation (420) are
evident in the recovered wavelet shown in Figure 28. The recovered
Example
Hamming Windowed Data Record
Figures 22 through 28
This group of figures illustrates the computation of the cepstra,
and wavelet recovery when the input data record is Hamming windowed.
Note the distortions present in the recovered wavelet.
The composite signal (256 points) is
y(nT) = x(nT) + .5x(nT30T)
where x(nT) = nT enT 0 n < 64
Figure Number Figure Title
22 Composite Signal
SNR = 40 dB
23 Unwrapped Phase Curve
24 Log Magnitude
25 Complex Cepstrum
26 Phase Cepstrum
27 Power Cepstrum
28 Recovered Wavelet
MSE = 4.24 x 103
65
o
I I 3 I 3 03 '
I I 3 I 1 6
I
:I 33
i' el
.II
:Ii
I lI I *I '
D: !
I 3 !
i i i *
1 1 11i
II I 1 .
II
J *J* *
! ,i 1 i
. 3 I a
SI i I I t I
* 3I I I 43
I I 3 I II I 4
' I 33 3 8
i 3 3 4 1 0,
! ! o Sn
I.~~ a a.
'3>!  ^ 3 ^L ^ .s 430
3 3 3 1 .4 4 M
* i ia .* ;? Z
* 3 .1 I ' *
SI* I I I I C
I1 I 4 I Ci
:::L
!i 1 .
I I I I I 
I I I I 3 *3 
I l i 3 a
3 3 I 3 I3
3 ~ ~ ~ I I. 3 '
I I I I
o o 3, *4, *3 K,
0 0
I I I
.J 4;r
3 3 3 3 4 .. I
66
o
a 1 a a'
e I I Ii I
e I l e !
I ! 
.i!
a a 1 a
a a ai
I I 4 I q
a a
a a a a
II tt i
I a i
i D
I I I o
a a a
!
i I i
1
%
. I I *
* .. C
*! I 1
I I a a a
a a a a
* a I a a WI
* a I I V
SII I
I a:
a a I a a Co
Soa
i i 1 
a I I a go
a* a a a i 
a a a a a a
a,
a a' a a. a a,
a I : a
I
I  I 0
'! *
< I 1 C
i* D
1 a a t
<^ ~ ~~ r ?" t
I II
4.. Ia, I 444 ta,
a a I I I 1*4 '
a~ I t
II atI :
I~ ~ a a aa
I a. I
a a I a ~
67
* 5 6 S I
:~ :
. I ,. 4 6 4

S* 4
I I I I t
i I 4 S S S
i S 44 6 S
I QU+ CCU~ 41UIII CIO~~ CI eLI
6 8 4 6 80
i .i I
I' i i
elI I
i . I I
SI 5 4' I 8
i i i I e
e i ii
I i i I
8 I 6 5
3 5 I I 4 5t
I 4. S 5
I I I I
SI S 8 I
I I I S I
S 4 1 5
5 8 4 4 5
* 3 4 I S
; *
Ie 6 8 I I
S 5 5 6 IS
I 8 46
I 5 4 6 I 5 C
*I 5 5 4
S4 S S S o
I 3 I 4 4 I 50
I 6 I 5 43
CI I,
6 6I
* 64 I SC
* I 81 50
I 5 5 4 I I .
* 5; S 4
* I *' 3 5I
0 I* 5 4 8^ 8
* I 5 4 6 6
* 8 485
* 4 S .
I
* aI a a I
! :
a a.
* .a *
* *b**.b~*a* t*ba S *
a a_ a a^ a. a a
a a a a a. a. a, a ^.r^ ^ m nii ^,,, l'L
a a
! a .
a a I I a. a
a a a a I
L
a a a a a.. aa
a a I I a a
a a a t
a aa
a a a a
a * a .
a..
a
a a a .. I
a a~ ac a a a aa
a~ a a a. a
a a a .
*,
a I a a .
a a
a a a a a
a a
a a a
:~ *.I
r*
: .
*
I
a 'a a' a a
a aI I a a
a a a a a. a
S* aaa
4 a a1 a. C aa
a aT a a a. a
a a a a 0 a0
a aI
I
69
: ".:
II I I .
* I I I
* S I I a0
* a I a__ c_ ^
* c
* S S I a c
* a oS
* S I S ,
c l $l r wr .usrrr
. S a a,
: i ie
aI I S0
a a a a a c
a a* aa cc c
* S I c
!
:
a a SN
S I S Sc
a a c,
* a a a
S :*
I I I Sc
* S I S 0 S
I Ia I
a. a I i
at.^titi^^._^ ^ ai.^t,n a^"** ** cap
a aI I 5 A
a. a a a
aI I c 0I
a a a a *O a0
<* a at a* a o
a a a a a c
a
70
a ... a~ u~~rr~ Il~~rrr
So
I
I S SI
S S
I I S a 0
; '"ro
I I S a a t n
;
1 :g
S ft
i
!
*_ S_ J S I_ I
II .
S I I 5
S1 o
* S S S I ft
0I S 0 S I
* I S t 5
* I S1 I S t
* I I 5 ft
* a I S f
* I I 5 f
SI I S a0
I SI I 5 ft C
* S I I S ft
*I I afp
* aS I f
* I S 5f
I S I 5
* I 55f
SI I I0
* I I 5 5 a
* I I 00
~ ~~ I 1
*~ S S If
* I S I 'I 0
* ~ ~ I I0
I S II 0: C
*I S S f
,c~c~ch ec,,,, r,,* I 5 I 5 ftc~e~w
* a I If
I I I~ I 0
* I I 5 f
I I I I
. o.
'oc.:
j i . i i'
a aI I **
a I SI
I S *e
I t<
8 S *! S
*,*,rt.>,,li , IIn*i Tw . u M ii
I S *
II SI *
_* a S I **
*I I *
P" ^r i 
*I I S **
a I
Ii S
*I a I C
i
~ ~ S I S
a S I
S S S a *lc~~~l r ~ccl c~~
S
a
*II S a *C
I a a I IC I
I I 8 C
I
S I .
aaC I C C 5 I I
I *C5a
a a
a Ia
I *
a aa
I
a a C a
C^ T* C "
0 0 0 D 0 0
,L~IUc~c~~,
cc~c~,U.
..
..
.
wavelet of Figure 28 is a result of smoothing only the peaks in the
complex cepstrum. When not only the peaks but several adjacent points
were smoothed the distortions clearly diminished. It is to be noted
from equation (418) that if E2 vanishes then the echo may be
recovered exactly;thus E2 is indeed an error term.
If B is less than one, a parallel analysis to that given above
predicts peaks at positive quefrencies, and reveals that a distorted
version of the basic wavelet rather than the echo will be recovered.
The theoretical analysis presented above gives considerable
insight into the wavelet recovery process when the window is
approximately constant over the duration of the basic wavelet.
A second analysis is presented below which gives further insight
into the effects of windowing on epoch detection and wavelet recovery
under different conditions.
Consider equation (49), it may be rewritten as
n Ll
Y(z) = X (z)+az o0 x(nTk T)w(nT)zn+E+E3 (421)
n=O
n L1
where E3 = az 0 x(nTk T)(w(nT+n T)w(nT))zn (422)
n=O
n
thus Y(z) = X (z)(1+az 0)+E+E3 (423)
Both El and E3 are considered error terms in that if they are not
present, then one could obtain the complex cepstrum, remove the
contributions due to the echo impulse train, perform the inverse
complex cepstrum operation, correct for the windowing by multiplying
the recovered sequence by (w(nT))1 and obtain the basic wavelet.
The error term El has already been discussed in the section on
the exponential window, and will again be assumed to vanish. A
close examination of the error E3 reveals it arises from two distinct
sources. E3 can be written as
no Lno1 n Ll
E3 = az no x(nTk T)[w(nT+n T)w(nT)]z +az o0 x(nTk T)w(nT)zn.
Sn= n=Ln
The first term in the equation above represents an error introduced
by the window because the window weights the basic wavelet and echo
differently, the second term is again the.truncation error due to
the fact that the echo may extend beyond the record under consideration,
and as in the section on the exponential window it is neglected.
This leaves only the.windowing error to be considered; obviously E3
vanishes when w(nT)=w(nT+n T). This again suggests the rectangular
window for wavelet extraction. Other common windows for should be
acceptable, provided x(nTk T)(w(nT+n T)w(nT)) is small. For
most commonly used windows (Hamming, Hanning, etc.) this requires
n to be small (i.e., the delay between the basic wavelet and echo
must be short). Carrying the analysis further we may rewrite
equation (423) as
E n
Y(z) = (X(z) + 3 )(1+az ) (424)
1+az
thus
E n
Y(z) = log Y(z)= log( (z)+  ) + log (l+az ) (425)
From equation (425) we see that the side of the complex cepstrum
From equation (425) we see that the side of the complex cepstrum
I
on which the echo peaks occur will be determined by the echo amplitude
(a) just as in the rectangularly windowed case.
Assuming that the term log (l+azn ) can be filtered out, we
can extract the wavelet by (assuming a
E
YRW(nT) = 1Z (X (z) +  )
l+az
= Z1(Xw(z) + E3(1azno+a2z 2n...)) (426)
YRW(nT) = w(nT)x(nTk T)+[ax(nTk Tn T)(w(nT+n T)w(nT))]
RW 0 0 0 0
*[6(nT)a6(nTn T)+a26(nT2n T)...] (427)
Multiplying (427) by (w(nT))l, we obtain
w(nT+n T) 2
YR(nT) =x(nTkT) + [ax(nT+k Tn T)( w(nT) )]a x(nT+k T2n T)
R o o o w(nT) o o
w(nTn T)
(1 w + .. (428)
Thus the basic wavelet is recovered though it is distorted by the
windowing error terms. A similar derivation applies if a>l, yielding
recovery of the echo rather than the basic wavelet, but with similar
distortion terms.
From the two analyses considered it is evident that there are
two cases when we may expect satisfactory wavelet recovery after
windowing with one of the windows commonly used to reduce leakage:
(1) When w(nT) is relatively constant over the duration of
the basic wavelet. This requires the signal duration to be short
compared to the total record length. The experimental results
presented for the signal of the form given in equations (44) and (45)
(with b=1.0) are an example of this case. Even though the actual
signal duration is a considerable portion of the total record length,
90% of the signal energy lies in a region only 3 or 4% of the total
record length. Thus the effective duration is quite short. The
first theoretical analysis presented accurately predicts the results
of these experiments.
(2) When the window w(nT) is approximately constant over the
w(nT+noT)
echo delay n that is w(nT 1. This requires n to be short
compared to the total record length. Additional computer computations
were conducted to verify the findings of the second analysis.
Utilizing the signal given in equations (44) and (45) with b=.1,
it is found that echo epoch detection and successful wavelet recovery
(in the no noise case) are possible for echo delays as short as 5
sample times when the rectangular window is used. When the same
signal is Hamming windowed,echo detection and satisfactory wavelet
recovery arepossible for small echo delays (n <10) though even for
n equal 5 the MSE is considerably greater for this window than for
the rectangular window.
It should be noted that no assumptions (other than regarding
the magnitudes of a.and B) were made in the analyses presented until
the filtering of the peaks in the complex cepstrum was considered.
Essentially both analyses are correct up to this point though the
nature and magnitude of the error terms differ. This leads to the
observation that when a1, and neither approximation is valid,
peaks should be present at both positive and negative quefrencies.
One set is predicted by the first theoretical analysis and the other
set is predicted by the second analysis. This is precisely what
was observed in many of the experimental outputs.
Conclusions
The exponential window
Windowing the input data with the exponential window was expected
from theoretical considerations to perform as well as the rectangular
window in the no noise case (and somewhat better if echo truncation
is a problem). Over the duration of the signal this is verified
by the experimental results at high SNR; however, correction for
the exponential window introduces some distortion into the recovered
record outside the signal duration. At low SNR the exponential window
performed better than the rectangular window. This is undoubtedly
due .to the fact that the composite signalexamined occurred at the
beginning of the data record and thus the exponential window tends
to reduce the noise content of the total record much more than the
signal content. Correspondingly when the composite signal occurs
near the end of the data record, the exponential window degrades the
recovered wavelet. It is similarly noted that while the echo epoch
peak in the complex cepstrum is diminished,.the detection threshold
is unchanged when the composite signal occurs near the beginning of
the record but increases when the composite signal occurs later in
the record. Finally it was noted that the exponential window may be
used to reduce aliasing of the echo impulse train, and to prevent
ambiguous peaks in the multiple echo case. Unless echo truncation
or one of the problems cited above is present, an exponential windowing
of the input data sequence appears to offer no advantages over the
rectangular window.
The common windows
Application of the Hamming, Hanning, and Tapering windows is
very detrimental to wavelet recovery. Satisfactory recovery is
not possible unless:
(1) the echo delay is small compared with the total record, or
(2) the duration of the basic wavelet is short compared with
the record length.
Even in these cases wavelet recovery is degraded. It is interesting
to note that in spite of the distortion caused by windowing, echo
epoch detection is always possible (at least in the no noise case)
and that while peaks maybe introduced on both sides of the complex
cepstrum, no shifting of the peaks from the proper echo epoch is
present. Windowing does however, increase the SNR at which echo
epoch determination is possible. This indicates that cepstral
techniques may be useful in resolving signals of similar, but not
identical, shapes. The theoretical analyses presented predict the
errors produced in the recovered wavelet by windowing with the
Hamming window quite well. A similar analysis may give some insight
into the distortion introduced in the recovered wavelet when the
echo is truncated.
WINDOWING THE LOG SPECTRUM
Since the complex cepstrum contains well defined peaks, it
seems reasonable that windowing the log spectrum be considered to
reduce the leakage due to these peaks. Only the effect of the Hamming
window on the log spectrum is studied, as this window seems repre
sentative of the windows used to reduce leakage and it is nonzero
at its endpoints making correction less difficult. After windowing
the log spectrum, the complex cepstrum is computed and the peaks
present due to the echo are removed ;this sequence is then inverse
transformed to obtain the recovered log spectrum which is (unless
otherwise noted) corrected by multiplying the recovered log spectrum
by the inverse of the windowing sequence.
The following observations are made on the effects of windowing
the log spectrum.
(1) A spreading of the peaks in the complex cepstrum is noted;
e.g., a single positive peak is reduced in amplitude and has two
adjacent negative peaks.
(2) Empirically, it is found that the "filtering" of the complex
cepstrum is critical. Smoothing the complex cepstrum at single peak
points as used previously proves to be totally inadequate, as the
recovered wavelet is virtually unrecognizable even when no noise is
present. When both the peak and the two points adjacent to it are
smoothed, recovery is possible. Figure 29 compares the MSE of the
recovered wavelet (when adjacent points are smoothed) to that
obtained when only rectangular windowing of the log spectrum is used.
As can be seen for the no noise case, windowing the log spectrum
results in about the same MSE as for the rectangularly windowed case,
but performance rapidly deteriorates when noise is added. The wavelet
becomes unrecoverable at a SNR of 14 dB.. The echo detectability
threshold (20 dB) is also increased by 12 dB over the rectangularly
windowed case.
* Rectangular Window
* Hamming Windowed Log Spectrum
MSE
106
10
7,
1 4
107 I I I I  4
10 20 30 40
SNR(dB)
Figure 29 MSE of the Recovered Wavelet when the Log
Is Hamming Windowed.
Spectrum
(3) Several runs have been made with no correction of the
recovered log spectrum, and it was found that wavelet recovery is
impossible if correction is not made.
Interpretation of the Results
Let us consider the effects of windowing the log spectrum in
the single echo case.
n
Y(z) = log Y(z) = log X(z)+ log (l+az ) (429.)
Windowing ?(z), we obtain
n
W(z)Y(z) = W(z) log (I+az ) + W(z) log X(z)
where W(z) is the window.
Assuming a
power series, and inverse ztransforming equation (429), we obtain
the complex cepstrum.
2 3
(nT) = w(nT)*(nT)nT)*(nT) T)*(a6(nno) 6(n2n )+ 6(n3n )...) .(430)
Thus, it is observed that windowing the log spectrum convolves the
complex cepstrum with the transform of the window. In the case of the
Hamming window this means the complex cepstrum is convolved with the
sequence
w(nT) = .226(nT+T)+.566(nT).226(nTT)
This accounts for the spreading of the peaks noted experimentally.
If the peaks due to the echo can be "filtered" from the complex cepstrum,
then we may obtain the basic wavelet by correcting the recovered
log spectrum and completing the remainder of the wavelet recovery
algorithm. It is observed that the contributions in the complex
cepstrum due to the echo are no longer isolated peaks but in fact
extend to the points immediately adjacent to the main peaks thus
we expect that "filtering" over a single peak will be inadequate
for wavelet recovery as is observed.
The rapid deterioration of the MSE with increasing SNR is prob
ably also connected with the filtering problem but this is not
well understood. The increase in the echo detectability threshold
is undoubtedly due to the fact that the peaks are reduced in height
by convolution with the sequence given above.
Conclusions
While windowing the log spectrum with the Hamming window reduces
leakage, it is apparent that other factors (principally the problem
with "filtering" the complex cepstrum) tend to degrade wavelet
recovery especially at low SNR. Thus windowing the log spectrum
with the Hamming window cannot be recommended as part of the general
wavelet recovery algorithm. The performance of other similar windows
is expected to be the same.
HANNING THE LOG SPECTRUM
It has been reported that the Hanning smoothing of the log
magnitude (that is, the real part of the log spectrum) results in a
decrease in both the MSE and the echo epoch detectability threshold
in the presence of additive noise [3]. Since this is closely related
to windowing the complex cepstrum, it is appropriate to discuss this
topic at this time. Smoothing consists of convolving the sequence to
be smoothed with a smoothing sequence (e.g., the Hanning smoothing
sequence is .25, .50, and .25). Smoothing is ordinarily used after
an FFT to reduce leakage, as this is equivalent to windowing the
input prior to the FFT with a Hanning window convolutionn in the
frequency domain is equivalent to multiplication in the time domain).
In the present context it serves a quite different function, smoothing
the log spectrum is essentially a frequency invariant low (short)
pass filtering operation. Consider the single additive echo case.
n
YS(z) = [log X(z) + log (1+az )]*W(2) (431)
where W(z) is the smoothing function or equivalently it may be
regarded as the impulse response of a filter. Inverse ztransforming
the above equation, we find
ys(nT) = ((nT)+e(nT))w(nT) (432)
where x(nT) is the complex cepstrum of the basic wavelet, e(nT) is
the train of 6 functions due to the presence of an echo, and w(nT)
is the inverse ztransform of the smoothing function.
For the Hanning weights given above, w(nT) = (l+cos ).
Thus we see that Hanning smoothing tends to reduce the high quefrency
NT
(near ) contributions to the complex cepstrum, and is equivalent
to windowing the complex cepstrum.
If only the log magnitude is smoothed, then
n n
Ys () = Re(log X(z)+ log(1+az O))*W(z)+jlm(log X(z)+ log(l+az o)
n n
= (loglX(z)l+ logll+az O)*W(z)+j(Phase X(z)+Phase (log(l+az 0))).
(433)
Evaluating the inverse ztransform on the unit circle z=ejmT and
noting that
ReY (z) = loglX(eJmT) + log(l+a2+2a cos amnT) (434)
is an even function of w, and
WT I a sin(n wT)
jImY(z) =jPhase X(ejT) +jatanl( cos(n T) (435)
is an odd function of w, we see that the transform of equation (434)
is the even part of the complex cepstrum y(nT) (i.e., it is
y (nT) = (y(nT) + y(nT)) and the transform of equation (435) is
the odd part of the complex cepstrum y(nT) (i.e., it is y (nT) =
I (9(nT)9(nT))). Thus the smoothed complex cepstrum is
9s(nT) = 9e(nT)w(nT) + 0o(nT) (436)
but
Se(nT) = n(9(nT)+y(nT))= (x(nT)+x(nT)+^(nT)+e(nT)) (437)
and
yo(nT) = ^((nT)9(nT)) =(x(nT)x(nT)+e(nT)e(nT)) (438)
where e(nT) is the train of impulses in the complex cepstrum due to
the echo. From the above analysis it is clear that the log magnitude
terms produce 6 functions on both sides of the origin, and the log
phase terms produce 5 functions which reinforce on one side of the
origin and cancel on the other to produce the onesided e(nT). When
the log magnitude is smoothed the 6 functions produced by it are
reduced by the window w(nT); thus they no longer cancel on one side
with the 6 functions due to the phase portion, and peaks may be
expected on both sides of the origin (though the peaks on one side
should be substantially greater). Thus we observe that a minimum
phase sequence (normally zero on the negative side of the complex
cepstrum) may produce contributions on both sides of the complex.
cepstrum if the log magnitude alone is smoothed. If both the log
magnitude and phase are smoothed their contributions will be reduced
equally and thus they will still cancel completely on the negative
side of the complex cepstrum.
Experimental Results
A Hanning smoothing performed on the log magnitude and phase
produces a windowing of the complex cepstrum as expected. For the
case considered (n =30) the echo epoch detection threshold appears
to be slightly lower than for the unsmoothed case. This is undoubtedly
due to the fact that the primary echo epoch peak is only slightly
reduced by the window while the peaks associated with noise located
at higher quefrencies are considerably reduced. Had the echo epoch
been at a higher quefrency (e.g., no=100) the detection threshold
would have increased greatly. As can be seen from Figure 30, the MSE
of the recovered wavelet is considerably reduced by the smoothing
process. Only in the no noise case did the unsmoothed process prove
to be superior.
Smoothing the log magnitude (only) produced similar results as seen
in Figure 30, but this is clearly inferior to smoothing both the log
magnitude and phase when additive noise is present. Since the power
85
Rectangular Window
Rectangularly Windowed Complex Cepstrum
10 a Hamming Smoothing of Log Magnitude Only
o Hanning Smoothing of Log Spectrum
103
MSE 104
0
MSE 10 0
105
10
107
107  I I I
10 20 30 40
SNR(dB)
Figure 30 MSE of Recovered Wavelet when .the Log Spectrum
Is Hanning Smoothed.
cepstrum is independent of phase information, the power cepstrum
in this case is identical to the power cepstrum produced above.
Kemerait [3] has reported observing echo peaks at both positive and
negative quefrencies only when a>l. Our results clearly contained
echo peaks at positive and negative quefrencies whether a>l or a
The positive peaks are dominant when a
dominant when a>l. This clearly confirms the theoretical analysis
presented. Since echo peaks are present at both positive and negative
quefrencies the complex cepstrum must be smoothed on both sides.
If only the dominant peaks are filtered the MSE of the recovered
wavelet is degraded slightly.
An attempt was made to improve the MSE performance by not
smoothing the log spectrum at low frequencies (where most of the
signal power is present) while continuing to smooth the higher frequency
components. This is again essentially a filtering process but in
this case the filtering is no longer frequency invariant. The results
from this type of smoothing process were found to be considerably
worse than those found above, and its use was discontinued.
Conclusions
The Hanning smoothing of the log spectrum is essentially a
frequency invariant low (short) pass filtering operation on the log
spectrum, or equivalently it may be regarded as a windowing of the
complex cepstrum. In this light the effect of smoothing on echo
epoch detection is easily understood, as it does in fact simply reduce
echo peaks occurring at high quefrencies more than those at low
quefrencies.
