Reference

NATL INST. OF STAND & TECH

NBS

Publications

AlllDb 3MD173

iiM- NBS TECHNICAL NOTE 1 047

U.S. DEPARTMENT OF COMMERCE / Notionol Bureau of Standards

Deconvolutlon of

Time Domain Waveforms

in the Presence of Noise

-QC1 100

U5753
NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards' was established by an act ot Congress on March 3, 1901. The Bureau's overall goal is to strengthen and advance the Nation's science and technology

and facilitate their effective application for public benefit. To this end, the Bureau conducts

research and provides: (1) a basis for the Nation's physical measurement system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety. The Bureau's technical work is per- formed by the National Measurement Laboratory, the National Engineering Laboratory, and

the Institute for Computer Sciences and Technology.

THE NATIONAL MEASUREMENT LABORATORY provides the national system of

physical and chemical and materials measurement; coordinates the system with measurement

systems of other nations and furnishes essential services leading to accurate and uniform

physical and chemical measurement throughout the Nation's scientific community, industry,

and commerce; conducts materials research leading to improved methods of measurement,

standards, and data on the properties of materials needed by industry, commerce, educational

institutions, and Government; provides advisory and research services to other Government

agencies; develops, produces, and distributes Standard Reference Materials; and provides

calibration services. The Laboratory consists of the following centers:

Absolute Physical Quantities' — Radiation Research — Thermodynamics and

Molecular Science — Analytical Chemistry — Materials Science. THE NATIONAL ENGINEERING LABORATORY provides technology and technical ser- vices to the public and private sectors to address national needs and to solve national problems; conducts research in engineering and applied science in support of these efforts; builds and maintains competence in the necessary disciplines required to carry out this research and technical service; develops engineering data and measurement capabilities;

provides engineering measurement traceability services; develops test methods and proposes

engineering standards and code changes; develops and proposes new engineering practices;

and develops and improves mechanisms to transfer results of its research to the ultimate user. The Laboratory consists of the following centers:

Applied Mathematics — Electronics and Electrical Engineering- — Mechanical

Engineering and Process Technology' — Building Technology — Fire Research —

Consumer Product Technology — Field Methods.

THE INSTITUTE FOR COMPUTER SCIENCES AND TECHNOLOGY conducts

research and provides scientific and technical services to aid Federal agencies in the selection, acquisition, application, and use of computer technology to improve effectiveness and

economy in Government operations in accordance with Public Law 89-306 (40 U.S.C. 759), relevant Executive Orders, and other directives; carries out this mission by managing the Federal Information Processing Standards Program, developing Federal ADP standards

guidelines, and managing Federal participation in ADP voluntary standardization activities; provides scientific and technological advisory services and assistance to Federal agencies; and

provides the technical foundation for computer-related policies of the Federal Government.

The Institute consists of the following centers:

Programming Science and Technology — Computer Systems Engineering.

'Headquarters and Laboratories at Gaithersburg, MD, unless otherwise noted; mailing address Washington, DC 20234. -Some divisions within the center are located at Boulder, CO 80303.
noMAL feasxAV or JiTAMDAXDS

UBKABT

OtC 1 1981

Deconvolution of ;;jr

"^ Time Domain Woveforms "TM

in ihe Presence of Noise ^^

N. S. Nahman

M. E. Guillaume

Electromagnetic Technology Division

National Engineering Laboratory

National Bureau of Standards

Boulder, Colorado 80303

M£.

U.S. DEPARTMENT OF COAAMERCE, Malcolm Doldrige, Secrerory

NATIONAL BUREAU OF STANDARDS, Ernesr Ambler, Direcror

Issued October 1981
NATIONAL DUP>EAU OF STANDARDS TECHNIGM NOTE 1 047

Nor. Dur, Srond. (U.5.), Tech. Nore 1047, 122 pages (Ocrober 1981)

CODEN: NDTNAE

US GOVERNMENT PRINTING OFFICE WASHINGTON: 1981
CONTENTS

Page FOREWORD V

1. INTRODUCTION - 1

2. GENERAL PRESENTATION OF THE DECONVOLUTION PROBLEM 4 2.1 An Exact Solution Exists for Exact Waveforms 4 2.2 Experimental Situation Yields Inexact Waveforms 4 2.3 The Regularization Operator 5 2.4 Possible Methods of Deconvolution 6 2.4.1 Model Fitting 6 2.4.2 Time Domain Processing 6 2.4.3 Transform Domain Processing , 6 2.5 Loss of Information in the Convolution Process 6

3. EQUIPMENT AND METHODS USED IN THE EXPERIMENTS H

3.1 Acquisition of Waveforms H

3.2 Interactive Processing 11

3.3 Qualitative Measurement of the Quality of One Solution 11 3.4 Quantitative Criteria 11

4. FREQUENCY (TRANSFORM) DOMAIN DECONVOLUTION 15 4.1 Introduction

4.2 Methods Considered 17 4.3 The Two Parameter Method 18 4.3.1 The Filter Structure 18 4.4 Application of the Two Parameter Method 19 4.4.1 Procedural Summary 19 4.4.2 The Estimation of the Insertion Impulse Response of 300 Meters of RG 58 C/U

Coaxial Cable 19 4.4.3 The Estimation of the Insertion Impulse Response of a Broadband Antenna ... 21 4.5 On the Nature of the Error, e(k) 22 4.5.1 Error in the Absence of the Filtering Operation and Computation Noise .... 22 4.5.2 Error in the Presence of the Filtering Operation and in the Absence of Computation Errors 22 4.5.3 Error in the Presence of the Filtering Operation and Computation Errors ... 25 4.6 The One Parameter Method 28 4.6.1 The Filter Structure 28 4.6.2 Conversion of a Step-Like Waveform into a Duration Limited One 35 4.6.3 Application of the One Parameter Method 37 4.6.3.1 Procedural Summary 37 4.6.3.2 Estimation of the Insertion Impulse Response of 300 Meters of RG 58 C/U Coaxial Cable 38 4.6.3.3 The Insertion Impulse Response of a Filter for a Step-like, 100 ps Transition Duration, Reference Waveform Generator 39
Page

5. TIME DOMAIN DECONVOLUTION 78 5.1 Introduction 78 5.2 Discrete Classical Time Domain Deconvolution 78 5.2.1 Convolution 78 5.2.2 Z-transform Representation 79 5.2.3 Deconvolution 80 5.2.4 The Effect of Errors in Classical Deconvolution 83 5.2.5 The Critical Dependence of Classical Deconvolution on the First Point in x(k). . 85 5.3 Discrete Iterative Time Domain Deconvolution 87 5.3.1 Iterative Method 87 5.3.2 The Effect of Errors in Iterative Deconvolution 89 5.4 Time Domain Regularization 90 5.4.1 Averaging Filter 90 5.4.2 Composite Filtering 91 5.4.2.1 Introduction 91 5.4.2.2 Filter Description 92 5.4.2.3 Frequency Domain Analysis of the Composite Filter 93 5.4.3 Application of Time Domain Filtering 94 5.4.3.1 Practical Procedure 94 5.4.3.2 Estimation of the Insertion Impulse Response of 300 Meters of RG 58 C/U Coaxial Cable 95

6. SUMMARY AND CONCLUSIONS 112 5.1 Summary 112 6.2 Conclusions 112

7. REFERENCES 114
FOREWORD

The work reported here is the result of a continuing collaboration between the Centre National

D'Etudes Des Telecommunications, France, and the National Bureau of Standards. From October 1978 through September 1979, Dr. Michele Guillaume was a Guest Scientist in the Time Domain Metrology Group of the Electromagnetic Technology Division, NBS , Boulder Laboratories, Boulder,

Colorado. At that time Drs. Guillaume and Nahman worked together to produce the experimental results

reported here. It has taken some time for the authors to complete their analyses of the experimental

data, to provide some analytic bases for understanding the deconvolution and filtering processes, and to prepare this manuscript. All of these latter activities were accomplished with the authors located in France and the USA, respectively. Thus, there has been some delay in publishing this work.

DECONVOLUTION OF TIME DOMAIN WAVEFORMS IN THE PRESENCE OF NOISE

N. S. Nahman and M. E. Guillaume t Electromagnetic Technology Division

National Bureau of Standards

Boulder, Colorado 80303

Deconvolution or inverse filtering was used to determine the impulse response

of a system using noisy input and output time domain sequences (discrete data)

.

Frequency and time domain methods were studied along with the synthesis of the filters required to obtain stable and smooth results. For the methods studied it was concluded that the superior technique was provided by an optimal frequency

domain method implemented via the FFT. Also, it is pointed out that the time

domain methods are only in their infancy and still retain the promise of avoiding

transform domain filtering. Examples are presented in which the impulse responses

are determined in the presence of varying degrees of noise for a coaxial transmission

line, a wave-shaping filter, and a broad-band antenna.

Key words: deconvolution; impulse response^ inverse filtering; noise contamination;

time domain; waveforms.

1. INTRODUCTION

The analogy between a physical linear system and the convolution equation is well known. If K(t)

is the input waveform to a linear device whose impulse response is h(t), the output waveform is described

using the convolution equation,

^(t) = fx(T) h(t-T)dr =J x(t-T)h(T) dx (1-1)

or in simplified notation

7(t) = x(t) * h(t) <;i-2)

where "*" denotes convolution. The Laplace transform of y(t) is given by

Y(s) = X(s) H(s) (^-3)

which shows that the convolution process transforms to a product operation in the s-domain. The deconvolution process is defined as the solution of the convolution integral equation (1-1) when either x(t)

or h(t) is unknown, i.e..

h(t) = y(t) (1/*) x(t) (1-^)

H(s) = Y(s) ^ X(s) (1-5)

t Centre de Microelectronique de Grenoble, Centre National d' Etudes de Telecommunications, Meylan, France.
x(t) = y(t) (1/*) h(t) (1-6)

X(s) = Y(s) T H(s) (1-7)

Note that (1/*) which denotes deconvolution transforms to a division operation. Both cases, i.e., the solution for h(t) or x(t) , are basically the same with respect to the deconvolution problem.

The above convolution integral relationship is defined for continuous waveforms. In the problems

considered here the waveform data are not continuous but are discrete; the data are obtained from a sampling process followed by an analog to digital conversion.

For the purposes of this present work, only discrete variables or sequences are considered. When

the question of error is considered, it is relative to the differences between one sequence and another.

The error due to representing a continuous function by a sequence is not considered here. However, it should be recognized that such an error (aliasing) can be very significant when spectral information is desired [1,2]

.

A sequence of N-values for x, x(0), x(l) , x(2) , ... x(N-l), is formally denoted as {x(k)} where the k-th number in the sequence is x(k)

,

X = {x(k)} (1-8)

For convenience, here the notation x(k) will denote the sequence of values for x. For discrete variables, i.e., sequences, the convolution and deconvolution equations (1-1) through

(1-7) are respectively replaced by

krl k-1

y(k) = 2. x(i)h(k-i) = Yl x(k-i)h(i)

i=o

(1-9)

i=o

y(k) = x(k)*h(k)

(^_-L0) Y(s) = X(s) H(s) ^_-^^^

h(k) = y(k) (1/*) x(k) (;L_3^2) H(s) = Y(s) V x(s) (1-13)

x(k) = y(k)(l/*)h(k) ^_^^)

X(s) = Y(s) ^ H(s)

(-L_.L5)

In the remainder of this report, for the most part, the analyses will be carried out in terms of discrete

variables. Where continuous variables are used, the distinction will be clearly made. However, in the figures, for convenience in drawing the discrete waveforms and spectra, they will be drawn as continuous

functions, but labeled with the appropriate discrete variables.

Numbers in brackets indicate literature references at the end of this paper.
In the experimental situation the imperfections of the acquisition device and the measurement

system noise which is added to the signal make the available data not representative of the original

continuous waveform but of something close to it. The convolution relationship still holds; but, in general, the deconvolution may no longer be stable. That is to say, the solution which comes from the direct brute-force application of the mathematical equations may not represent the desired signal.

Furthermore, an exact solution does not exist and must be replaced by an estimate. The problem is to find an approximate solution which is stable under small changes in the initial data. To do that one uses a regularizing operator or filter. This filter depends on the characteristics of the processed waveforms and on the method used to perform the deconvolution. Particular emphasis is given to the fact that the experimenter controls the operation and is the judge of the quality of the deconvolution result.

In this report after a general presentation of the deconvolution problem, two possible solution

methods are examined. In each method the reasons for the instability are inspected and lead to the synthesis of the corresponding regularizing function. The experimental results include the evaluation

of the quality of the waveform obtained as a solution of the deconvolution process. Specifically, this report is divided into six chapters, the first being the present introduction.

Chapter 2 is devoted to a general presentation of the deconvolution problem and possible approaches for solution.

Chapter 3 describes the experimental equipment used to obtain waveform data and the various experimental quantitative criteria which may be used to measure the quality of a solution to a deconvolution

problem.

Chapter 4 discusses frequency domain deconvolution methods and presents detailed examples of solutions to deconvolution problems using frequency domain processing which include the synthesis of regularization filters.

Chapter 5 discusses time domain deconvolution methods and presents several examples, using time

domain processing, one illustrating a time domain solution to a deconvolution problem which was also

solved in Chapter 4 using frequency domain methods. The synthesis of regularizing filters for time

domain deconvolution is also discussed.

Chapter 6 presents a summary of the deconvolution methods considered and some general conclusions

about the solution to the deconvolution problem.
2. GENERAL PRESENTATION OF THE DECONVOLUTION PROBLEM

2.1 An Exact Solution Exists for Exact Waveforms

Consider the linear system shown in figure 2-1; if the input waveform x (k) and the measurement

system impulse response h (k) are assumed to be exactly known, then the exact output waveform y (k) is also known and is given by the convolution equation,

y^(k) = x^(k) * h^(k). (2.1-1)

Furthermore, if h (k) is the unknown, the deconvolution

h^(k) = y^(k) (1/*) x^(k) (2.1-2)

is a unique solution which exists only because x (k) and y (k) represent exact, error free waveforms

(sequences)

.

2.2 Experimental Situation Yields Inexact Waveforms

Practically, the exact waveforms x (k) and y„(k) are not available from measured experimental data,

because the exact values are corrupted by system noise and sampling errors. Consequently, the experimental data x(k) and y(k) replace the exact data x (k) and y (k) respectively. However, each representation, x(k) or y(k), belongs to a family of waveforms defined as "close" to the exact waveform; that is to say, x(k) belongs to the family [X],

X £ [X], (2.2-1)

such that for all elements x of this family,

P(x, X ) < 6 (2.2-2)

o — X

where p being defined for any kind of waveforms w and w by the relation.

'"1 ^ 21 P(w.%) = j^N ^ ^"i - "oi^

1/2' (2.2-3)

The parameter 6 may be interpreted as a measure of the error which is made when the waveform w is represented by w. Then for

X £ [X], p(x,x^) < 6 (2.2-4)

y £ [Y], p(y,y^) < 6^ (2.2-5)

These two relationships define families of waveforms such that "small changes" are allowed between the elements of each family.

The direct use of the deconvolution equation (1-3) when applied to the y(k) and x(k) waveforms

leads to a result h(k) which, in general, is not close to h (k) the exact solution, i.e., if [H] is the family of waveforms close to h (k) , then

h(k) i [H] (2.2-6)

and

p(h, h^)> 6^ (2.2-7)

hence, the solution is not stable under small changes in the initial data.
A strict mathematical solution is just not possible because, in general, the result is not stable [3] The problem is then to find an approximate solution stable under small changes in the initial data. The

requirement of stability has to be satisfied since it is associated with the physical determinancy of the phenomenon described by the convolution equation. A rule has to be chosen for the definition of the family [H] of possible solutions. In practice this choice is arbitrary; consequently, one needs to use supplementary information concernin the solution.

However, it is known that the convolution itself is stable under small changes in the initial data.

If h(k) is the obtained approximate solution, its convolution with the known input waveform x(k) must

give a result which belongs to the family [Y] of waveforms close to y (k) . If

c(k) = h(k) * x(k) (2.2-8)

and

p(c,y) < 6 ; c £ [Y] (2.2-9)

then

e(k) = y(k) - c(k) (2.2-10)

defines an error criterion whose properties will be used to judge the quality of the representation of h (k) by h(k). Moreover, it may be assumed that the approximate solution has to be smooth because the solution represents a continuous deterministic signal or system (in sampled data discrete form.) This

qualitative information will be used in addition to the error criterion to control the process which

leads to the accepted solution.

2.3 The Regularization Operator

Another way which employs linear algebra, i.e., matrix methods, may be used to describe the convolution process. The convolution

y(k) = h(k) * x(k) may be written as y(k) = A x(k) (2.3-1)

where A is an operator which maps the family [X] into the family of waveforms [Y] . It has already been

pointed out that the inverse operation corresponding to the operator A , which was expected to map [Y] in [X] , is not stable. Hence, A has to be replaced by a new operator R which really maps [Y] in [X] such that

x(k) = R y(k). (2.3-2)

R is called the "regularized operator". Its expression depends upon the properties of the two families of waveforms [X] and [Y] . These properties include the precision in the knowledge of a waveform when compared to the exact waveform (parameter 6) and the shape of the waveform itself. On the other hand, the method used to perform the deconvolution determines the practical effect of the noise

and influences the definition of the regularized operator. From the user's point of view, the properties of R must be easy to modify in such a way that a progressive improvement of the regularization

finally leads to the accepted estimate solution.

Functionally,

R = R {[X], [Y], a(6)} (2.3-3)

where a is a set of parameters which are used to take into account the different possible noise figures.
2.4 Possible Methods of Deconvolution

There are many different methods for solving the deconvolution equation [4,5]. In this section

three methods are briefly discussed in order to provide a broader perspective of deconvolution methods.

In the main part of this report only two of the three methods will be discussed. In practice, all

methods of deconvolution (including the three) have a common feature in that the amount of information

that can be extracted from the empirical curves (physical data) is limited. Consequently, some form of approximation is a necessary part of the solution.

2.4.1 Model Fitting

Model fitting consists in finding an analytical expression to represent the exact waveforms y (k) and X (k) before the deconvolution is performed e.g., Prony's Method [6,7]. In this case, from the inexact available data it is necessary to obtain a mathematical expression which is assumed to represent

the exact solution. This operation has to be made for each waveform with the additional condition that

the two obtained waveforms have to be compatible in the convolution process, i.e., that y' (assumed exact representation of y(k) must be the image of x' (assumed exact representation of x(k)) in the convolution operation. This method, model fitting, is not examined in this present experimental study.

2.4.2 Time Domain Processing

The problem is solved using only the time domain representation. The regularizing operator is also

defined in the time domain. This class of methods includes the inversion of the convolution matrix and

the iteration technique.

2.4.3 Transform Domain Processing

Transform domain processing methods include the method of moments (statistical) [4], the homomorphic

deconvolution [8] and the Fourier transform deconvolution. The latter method is considered here. Generally, in transform domain processing the available data are transformed into another representation

(moments, frequency, etc.) and the regularizing operator is defined in this new domain before the return to the time domain.

2.5 Loss of Information in the Convolution Process

The details of the observed waveform are smoothed by the convolution with the apparatus function.

It is impossible to distinguish between the effects of the noise and very small irregularities due to the signal itself. Moreover, the amount of information available in the output waveform must be com- parable to the amount of information needed for the description of the input waveform (or impulse

response) . This last point is very important and may be explained using a frequency domain presentation.

First of all, when the duration of the unknown impulse response waveform is smaller than that of the input waveform, figure 2.2, the corresponding discrete frequency domain spectrum magnitudes are as shown

in figure 2.3. It is evident from figure 2.3 that the discrete spectra X(n) and Y(n) do not contain

enough high frequency components to allow a complete description of H(n). It is then impossible to determine H(n) or h(t) unless additional information is used.
T'Jhen the duration of the unknown impulse response is larger than that of the input pulse, x(t) and

y(t) contain enough information to determine h(t), figures 2-4 and 2-5. Consequently, a possible waj' to remove the fundamental difficulty illustrated by figures 2-2 and 2-3, is to use an input signal x(t) of lesser duration. When the input pulse duration is less than that of the system impulse response, the only problem which remains is to find a way to manage the effects of the noise which always is present.

Such management is called a regularizing process or regularization.

In conclusion of this chapter on the general presentation of the deconvolution problem, it is important to reiterate (in slightly different words) a point made in an earlier study [5]

"A strict mathematical solution cannot be obtained

because of the ever present noise in the original

data, but the correlation between the mathematics

and the physical data must alwa^'s be kept in mind

to reach an acceptable estimated solution".

The effect on the final solution of the uncertainty of the initial data due to noise will depend upon the method used for the deconvolution. Two methods will now be studied in detail: the Fourier

transform processing method and the time domain processing method. In the frequency domain, it will be seen that the disturbing effect of the noise appears in the high frequency range. In the second method,

i.e., in the time domain method, it will be seen that the noise increases the value of the difference

between two consecutive samples, and this random deviation from the exact values causes instability in

the solution. Furthermore, the time domain method is very sensitive to the D.C. value or offset of the waveform; this leads to difficult^' in using the method.
X (k) oInput

Linear System

h (k)

o Output

Figure 2.1 Linear System Possessing an Impulse Response h^(k).
x(k) A h(k)

>k V

y(k)

N-1

Input

N-1

System Impulse Response

>. k

N-1

Output

Figure 2.2 The unknown h(k) convolved with x(k) where the duration of

h(k) is less than that of x(k).

Figure 2.3 The frequency domain spectrum magnitudes corresponding

to x(k), h(k) and Y(k) of figure 2.2.
x(k) h(k) y(k)

Input

—1-> k Kl ZS*. ^

k

N-1 N-1

System Impulse Response

N-1

Output

Figure 2.4 The unknown h(k) convolved with x(k) where the duration

of h(k) is greater than that of x(k).

N/2 N/2 N/2

Figure 2.5 The frequency domain spectrum magnitudes corresponding to

x(k), h(k) and y(k) of figure 2.4.

10
3. EQUIPMENT AND METHODS USED IN THE EXPERIMENTS

3.1 Acquisition of Waveforms

A minicomputer controlled sampling oscilloscope system was used to implement the deconvolution

experiments. The specific system employed was the NBS Automatic Pulse Measurement System (APMS) [9,10]. The oscilloscope bandwidth is very large (dc-18 GHz) and permits the observation of picosecond domain

signals. A block diagram of the APMS is shown in figure 3.1.

In operation, a repetitive signal synchronized to the oscilloscope time base is applied to the sampling head. The time at which the samples are taken is controlled by the minicomputer (using the digital to analog converter). The amplitude of each sample is digitized and stored in the memory. The effect of the noise which always accompanies the signal may be reduced by using either analog or digital

averaging techniques. The APMS can only acquire a limited finite amount of data (points); consequently,

a time window of finite extent is used to observe the waveform. The acquisition process is fast and provides a large number of data points (102A for example) which are stored in an array whose format is well adapted for further digital processing. The CRT terminal provides a powerful tool which allows the

operator to implement, control, and monitor signal acquisition, conditioning and deconvolution processing.

3.2 Interactive Processing

Because no exact solution exists, it is difficult to make the deconvolution process entirely

automatic. For that, one would have to use a general mathematical criterion whose definitions could

hardly take into account each particular configuration of available data. The resulting solution would

not be so precise as the one which is obtained if the judgment of the operator is taken as a part of the process. With the operator involved in the process the graphical capability of the CRT display is extensively used to maintain the contact between the operator and the digital processing. Each time a decision is made its effect appears in the display permitting a measurement of its effectiveness or suitability.

3.3 Qualitative Measurement of the Quality of One Solution

The operator uses all the available qualitative information. He may have a preconceived idea of the shape of the expected waveform. In this case he may recognize that the obtained result is close or far off the solution he expects. Often, the only usable criterion is the smoothness of the deconvolved

waveform which results from the deconvolution process.

3.4 Quantitative Criteria

Because different methods may be used for the achievement of the deconvolution, the quality of the result has to be measured using a quantitative error function. This criterion permits a choice between

several possible solutions. The error function [e] is defined by the relationship, figure 3-2.

[e] = [y] - {[d] * [x]} (3.4-1)

[y] and [x] are the two known waveforms (output and input), and [d] is the approximate solution of the deconvolution

[d] = [y] (1/*) [x]. (3.4-2)

11
If the deconvolution process is perfect (exact solution) all samples of the error function [e] will be

equal to zero. The use of this error function is justified by the following remarks. Remember that the convolution

operator is stable under small changes in the initial data. If the obtained result [d] belongs to the family [H] of waveforms close to the exact solution [h ], the result [y'] of its convolution by the input

waveform must belong to the family [Y] of waveforms close to the exact output. But the exact output [y ] is not known. Then [Y] is used as the reference for the definition of the error function (3.4-1). The error function must be carefully interpreted. Four parameters will be used for its characterization, the mean value e, the standard deviation a, the maximum value e and the minimum value e . . max mxn

The mean value is given by

N1

k=l J E e(k) (3.4-3)

Ideally e should be zero, hence in practice a value as close to zero as possible is sought. The standard

deviation is given by

N _ - E [e(k) - e]' k=l1/2

(3.4-4)

and must be small when compared to the maximum value of the reference waveform [y] . A typical acceptable

value for a is when it is the same order of magnitude as the noise on the output waveform [y] . Experimentally, it has been shown that the maximum and minimum values of the error, e and e . occur for max mm

some samples associated with the largest values of [y] , figure 3.2; also, for a well chosen regularizing

function le and e . will be of the same order of magnitude. Finally, it should be kept in mind

' max '

' mm

'

_

that in the course of the deconvolution process, f ' the goal 1= is to make the error parameters f e, a, e and max

e . very small compared to the maximum value of the reference waveform [yl. mxn

12
FLEXIBLE

DISK DRIVE

[Tt7|—

I PTPI

SYSTEM

MINICOMPUTER

HARD COPY

UNIT

VOLTAGE AXIS MEASURED SIGNAL

-•—I 14 BIT A/T]-*-

DIGITAL L

INTERFACE I

U BIT D/A

CRT GRAPHICS

TERMINAL

TIME AXIS SHEEP SIGNALH

WIDEBAND

SAXPLIKG OSCILLOSCOPE

SAMPLING

HEAD

Figure 3.1 The NBS Automatic Pulse Measurement System

(PTP and PTR are paper-tape-punch and

paper-tape-reader, respectively)

.

13
d(k)*x(k)

max

e(k)

mm

k

TX^Ay<V»*. .^. /V>»wii r-n. mmA ^ m ! ^ k

Figure 3.2 The error e(k) is the difference between the reference

function y(k) and the estimated function d(k) * x(k).

14
4. FREQUENCY (TRANSFORM) DOMAIN DECONVOLUTION

A.l Introduction

Frequency domain deconvolution methods operate indirectly on the time domain sequence x(k) ; they

operate on the transform sequence X(n) which is the Discrete Fourier Transformation (DFT) [11,12] of the time domain sequence. The DFT and its inverse IDFT are defined in the following way. The sequence X(k),

i.e. ,

x(kT) = x(0), x(T), x(2T), . . . x[(N-l)T] (A. 1-1)

is sequence of values spaced apart by an interval T. The DFT of x(k) is a sequence of complex values

given by

N-1 X(n Q) = Yl ^(kT) e~^^ '^''", n = 0, 1, 2 ... N-1 (A. 1-2)

k=0

where each value is spaced apart by an interval Q, where fiT = 2'n'/N. The inverse transformation IDFT is given by

N-1

x(kT) = ^ Z) X(n Q) e^ ^

'^^'^

, k = 0, 1, 2 ...N-1 (4.1-3) n=0

As indicated above, the notation x(k) and X(n) implies x(kT) and X(n Q) , respectively.

Errors introduced into the deconvolution process stem from three sources:

1. In the time domain - erroneous or inconsistent data and computation

errors in computing the error.

2. Transform Domain - filtering and computation errors.

3. Transformation Operation - computation errors.

As summarized in Chapter 2, noise or errors in the time domain data cause the data to be inconsistent in the sense that the input x(k) and output y(k) data (i.e. the two known quantities) are not re- lated to each other through a causal stable impulse response. For example, if for a given causal stable

system x (k) , y (k) , and h (k) are the input, output and system impulse response sequences, respectively,

then the sequences are uniquely related by the convolution equation (1-8)

y (k) = x (k) * h (k) (4.1-4)

o o o

Consequently, x (k) and y (k) together provide the data to exactly characterize the system h (k)

figure 4.1. Now assume that in the measurement of x (k) and y (k) the independent noise sources, n and n , corrupt x (k) and y (k) , respectively, to yield x(k) and y(k), figure 4.2. Clearly, x(k) and y(k) are different from x (k) and y (k) , respectively. Furthermore, in general x(k) and y(k) are not related

to each other through (4-4) and cannot provide h (k) through the deconvolution process as specified by

h (k) = y (k)(l/*) X (k) (4.1-5)

o o o

which specifically employs x (k) and y (k) and not x(k) and y(k) , respectively. However, the impulse

o o response

h(k) = y(k)(l/==--) x(k) (4.1-6)

can be calculated. This result may be stable or unstable, and if stable, it may resemble h (k) or a very noisy h (k)

15
The impulse response of a physical system is causal (zero for time less than zero, the time at which

the impulse is applied) and stable (returns to zero after the passage of time), d (t) in figure 4.3. By operating on the data sequences with the DFX the discrete Fourier transforms of x(k) and y(k) can be obtained; these transforms are the complex sequences,

DFT[x(k)] = X(n) (4.1-7)

DFT[y(k)] = Y(n) (4.1-8)

Also, it is true that the exact processes (4.1-4) and (4.1-5) have the DF transforms

Y (n) = X (n) H (n) (4.1-9)

o o o

H (n) = Y (n)/X (n) (4.1-10)

o o o but the ratio Y(n)/X(n) still does not represent H (n) , i.e.

o H (n) i Y(n)/X(n) (4.1-11)

Furthermore,

which leads to the conclusion

H(n) = Y(n)/X(n) (4.1-12)

H(n) i H (n) (4.1-13)

o Consequently, the DF transformation in itself does not remove the errors that change X (n) and Y (n) into X(n) and Y(n), respectively.

The amount of error in erroneous or noise corrupted data depends upon the signal to noise voltage

ratio, SNR, e.g., [x(k)] /a , where O is the standard deviation (3.4-4) of the zero-mean measurement max X X system noise present while x(k) was being measured. For a very large SNR, the noise has a very small

effect upon the resultant (sampled) data sequence. With decreasing SNR's, the noise becomes increasingly

dominant. Typically, for the most part, the noise spectrum is uniformly distributed (constant spectrum

amplitude) and is weighted or filtered by the network properties of the sampled-data measurement system. The net result is that the higher frequency components of a low pass signal well within the measurement

system passband are increasingly corrupted by noise with increasing frequency. Consequently, the corrupted DFT's X(n) and Y(n), respectively, with increasing corruption at the higher frequencies.

Furthermore, their ratio gives the corrupted transfer fuction H(n) , (4.1-12), whose IDFT yields the corrupted impulse response Tymi

I X(n) I h(k) = IDFT [H(n)] = IDFT^^ (4.1-14)

I- J

which may or may not be stable, and which may not even resemble the exact or true impulse response sequence h (k) : stable and unstable responses are shown in figure 4.3.

o

-J To reduce the influence of noise at the higher frequencies on H(n) and ultimately on h(k), H(n) is filtered with a low pass filter, R(n), which is judiciously chosen so as to provide a stable and smooth

impulse response, d(k). R(n) operates on H(n) to give D(n),

D(n) = R(n) : H(n) (4.1-15)

where : denotes that R(n) operates on H(n) . In the time domain,

d(k) = IDFT D(n) (4.1-16)

The error e(k) due to the estimate d(k) is then obtained from (3.4-1) and characterized as described in

Section 3.4.

16
Too much filtering will yield an impulse response which varies too slowly, while too little

filtering may yield an unstable result, figure 4.3. The choice of the low pass filter parameters is in fact arbitrary because the filtered result giving the impulse response, d(k), is subject to the interpretation of the viewer. Accordingly, the filtering process introduces an intrinsic error even though

filtering may be necessary to achieve a stable result, i.e., a causal stable impulse response.

In addition to the noise corruption of the time domain sequence and the transform domain filtering

distortion, the DFT , IDFT, and other computations can introduce numerical errors which also affect the choice of the estimate d(k). The DFT and IDFT are usually implemented using one of the fast Fourier

transform (FFT) algorithms [13]. DFT, FFT and other operations can be analyzed to estimate the resultant

signal to noise ratios where the noise is due to numerical truncation and rounding required by finite

computing register length [14]. However, the errors due to computations are relatively small compared

to those due to noise and filtering.

In summary, the overall deconvolution processes for time domain sequences using transform domain

filtering possesses three sources of error whose effects can only be reduced but not eliminated. Hence,

it is clear that the deconvolution process is inexact, i.e., the true answer, h (k) cannot be obtained

as has been pointed out earlier (Section 2.5), and reiterated here. However, because "hope springs eternal in the human spirit" techniques have been and continue to be developed by which x(k) and y(k) are operated upon in order to "determine" h (k) . The job cannot be done in an exact sense, that is, h (k) cannot be determined, but it can be estimated; the result is an

o estimate for h (k) . In many practical situations the results are good enough to be very useful. To be good enough and useful, the results reasonably approximate some preconceived view of the system behavior

derived from a physically-based mathematical model. The argument can be made to appear more plausible

by saying that the effects of the computation noise and the noise components n (k) and n (k) , figure 4.2

are relatively small; but a little noise can ultimately be just as devastating to the deconvolution

process as a large amount of noise. In fact, in practice, techniques for noise reduction or enhancement

of the signal to noise ratio are first used; then, that noise which remains will still be potent enough

to prevent exact deconvolution. The direct effect of noise can be seen in the analysis of time domain

deconvolution processes; analyses of the effect are presented in Chapter 5.

4.2 Methods Considered

Two frequency domain methods are considered in the remainder of this chapter; they are:

1. A two parameter filtering method

2. A one parameter filtering method.

The common objective of the two methods is to filter out the noise which inhibits the deconvolution

process from producing a smooth stable result. The filtering process is called a regularizing operation

and the filter, a regularizing filter, is denoted by R(n)

.

In the first method there are two parameters which fix the characteristics of the filter. Its simple mathematical form enables the viewer to easily envision and choose a trial set of values for the

two parameters based upon spectrum amplitude data in dB, 20 log|H(n)|, versus frequency n Q. After the

two parameters are selected, the viewer observes the time domain response (IDFT), d (k) . of the filtered

H(n) , D(n). The choice of the filter structure in this method was not based upon any strong a priori

conditions, but rather upon the weak condition to attenuate the noisy portion of the spectrum for H(n)

and thus reduce the error (3.4-1).

In the second method there is only one filter parameter; its value governs the gradual on-set of filtering to provide a smooth time domain response. Again, the viewer selects a trial value for the filter parameter and observes the time domain response, d (k) of the filtered H(n) , D(n) . However, in

17
this method the filter structure is based upon the a priori conditions to minimize a performance measure consisting of a weighted sum of the error power and the smoothness power (power being the energy per period or time window.)

4.3 The Two Parameter Method

4.3.1 The Filter Structure

The filter specifically operates on the real and imaginary parts of H(n)

,

H(n) = Hj^(n) + j H^Cn) (4.3-1)

D(n) = Dj^(n) + J D^(n) (4.3-2)

to yield

As it was indicated for (4.1-14), the notation R(n) :H(n) denotes that R(n) operates on H(n) which in general implies that the operation need not be the product R(n)H(n) . In the two parameter method, the operation is not a product and operates on H (n) and H (n) to yield

K 1

D,(n)

H„(n), < n < n

K — o

|H(n )| exp(-Mn), n ^ n £ N-1 (4.3-3a)

H^(n) ; < n < n

I — o Dj(n) =

I (4. 3- 3b)

o , n < n < N-1

o — —

The magnitude of D(n) is given by

[H^(n) + \i\{n)]^'^ , o < n < n^

|D(n) j = { (4.3-4)

lH(n 'o')| exp(-Mn) "^ , n < n < N-1 o— —

and

20 log[H^(n) + Hj(n)]^''^ , < n < n^ D(n)|,^ = { (4.3-5)

' 20 log|H(n )1 -20(loge)Mn , £ n < N-1

I The constant M is the cut-off slope of filtered H(n) , i.e., [D(n)i^g, figure 4.4 and is given by

-100 + iH(n ) I -2_dB (4.3-6)

(A-l)n^

The rationale for the structure of the filter was to provide a simple means for attenuating

magnitude lH(n)| response for n greater than some value n^ (i.e., in the noisy region), while

simultaneously totally eliminating the noisy phase function, e(n),

. Tan"^[H fn)/H^(n)], < n < n

)(n) = ) ' ^ - o (^.3_7)

, n < n < N-1

o — —

18
4.4 Application of the Two Parameter Method

4.4.1 Procedural Summary

In practice, the viewer inspects the data H(n) and selects the acceptable range, n , of the data

do o (relatively noiseless data); this specifies n and |H(n )|. Then, the viewer selects the cut-off co- efficient, A, which then completely specifies the cut-off slope, M, and D(n). Finally, the viewer calls

for the IDFT of D(n) which provides the estimated or deconvolved impulse response, d(k). If the resultant d(k) and its associated error criteria, (3.4-1) through (3.4-4), are not satisfactory in the judgment of the viewer, the process is repeated using different values of n and A. The procedure starting with the time domain sequences x(k) and y(k) and ending with the deconvolved

time domain impulse sequence d(k) progresses as follows:

1. Offset correction for x(k) and y(k). Eliminate the offset, i.e., insure

that x(o) = x(N-l) = and y(o) = y(N-l) = 0. If x(k) and y(k) are steplike waveforms so that x(N-l) and y(N-l) are not zero, they must be transformed to impulsive waveforms. See subsection 4.6.2.

2. DF transform x(k) and y(k) to obtain X(n) and Y(n); display |x(n)| and

dB

l^(-)ldB3. Calculate the complex ratio H(n) = Y(n)/X(n); display |H(n)| . an

4. Select the regularization filter coefficients n and A; then filter

o H(n) to yield D(n), (4.3-3) display |D(n)L^.

dB

5. IDF transform D(n) to obtain the time domain deconvolved result,

d(k). (3.4-2), (4.1-16); display d(k).

6. Calculate the error criteria (3.4-1), (3.4-3), and (3.4-4) display e(k) and the values for e , e . , e, and a. max . mm

.

7. Judge the quality of the deconvolved impulse response estimate, d(k),

in terms of the error criteria; if unsatisfactory, return to step 4 and change the filter parameters.

8. If the quality is satisfactory, print the d(k) graphic display,

related error criteria, and the numerical values.

Examples of the two parameter deconvolution method are given below for the determination of the impulse response of (1) 300 meter coaxial transmission line and (2) a very broadband antenna.

4.4.2 The Estimation of the Insertion Impulse Response of 300 Meters

of RG 58 C/U Coaxial Cable

In this example the insertion input and output waveforms of a 300 meter length of RG 58 CU coaxial

transmission line, x(k) and y(k), respectively, were recorded, figure 4.5; the signal to noise ratios

for these waveforms were relatively large, SNR-x 75 dB and SNR-y 46 dB. Deconvolution to obtain the insertion impulse response was then performed using varying degrees of regularization with the two parameter

filter. The DFT of the insertion impulse response is the scattering insertion parameter S^_(n). [17].

The discrete transformation and inverse transformation operations were implemented using the FFT and IFFT algorithms.

Next, various combinations of SNR-x and SNR-y were created by adding noise (a pseudorandom

sequence) into the original high SNR data recorded above. For each combination of SNR's, two parameter

deconvolution was performed, again using varying degrees of regularization.

19
Specifically, the following deconvolution experiments were performed:

1. SNR-x: 75 dB, SNR-y : 46 dB

No regularization: (N-1) Q = 250 MHz Moderate regularization = n fi = 100 MHz, A = 5 Strongo regularization: o n Q = 74 MHz, A = 6 o

2. SNR-x: 75 dB, SNR-y: 36 dB

No regularization: (N-1) fi = 250 MHz Strong regularization: n fi = 74 MHz, A = 6

3. SNR-x: 45 dB, SNR-y: 46 dB

No regularization: (N-1) Q = 250 MHz Moderate regularization: n Q = 100 MHz, A = 5

4. SNR-x: 40 dB, SNR-y: 40 dB

No regularization: (N-1) Q = 250 MHz Strong regularization: n Q = 74 MHz, A = 6

The results are summarized in figures 4.6 through 4.8 and more detailed information is presented in figures 4.9 through 4.23. Figure 4.6 brings together the results for the no-regularization cases of each of the experiments.

In each of these cases, H(n) was directly inversely transformed by the IDFT to obtain the time domain sequence h(k). Also, shown in figure 4.6 are the |H(n)| corresponding to the various h(k) responses.

do Generally, it is seen that the noisier a |H(n)| is, the noisier is the corresponding h(k). FurtherdB more, a change in the SNR-y of the output waveform is more critical than that of the input waveform with

regards to the noise characteristics of |H(n)| and h(k). For example, a 10 dB reduction in SNR-y

do yields very noisy results, figure 4.6b, while a 30 dB reduction in SNR-x yields a much smaller noise

increase in |H(n)| and h(k), figure 4.5c. This follows from the fact that in the latter case the

dB transmission properties of H(n) attenuates the effect of the noise components introduced at the input.

Figure 4.7 brings together the results for the regularized cases of each of the experiments. In each of these cases, H(n) is operated on by the two parameter (n ,A) regularization operator to yield

D(n), and the result is inversely transformed by the IDFT to obtain the time domain sequence, d(k). The horizontal scale is expanded to 50 ns/Div. (as compared to 200 ns/Div. in figure 4.6) in order to show the fine structure of d(k). Also, h(k) for the largest SNR case (figure 4.5a) is shown in ex- panded time scale as figure 4.5a for comparison purposes.

It is evident that the regularization or filtering process produces much smoother estimates, d(k),

of the actual coaxial cable impulse response, h (k) , than do the various non-regularized h(k). The

smoothest non-regularized response is shown in figure 4.7a while the other non-regularized responses

of lesser smoothness are shown in figure 4.6. The effect of regularization is very dramatic in reducing the noise from the impulse response. Before commenting any further on figure 4.7, it is necessary to consider the input and output

waveforms of figures 4.9 and 4.10, respectively. First, from figure 4.9, it is clear that the input

pulse waveform abruptly rises, at say t^ , and is zero prior to t^; also, the pulse duration was of the order of 10 ns . Second, from figure 4.10, the output pulse waveform possesses some small variations

prior to the large abrupt rise of the pulse waveform; also, there is a small variation from the smooth

decay of the pulse waveform at about 800 ns. The main rise and subsequent return to zero of the output waveform occurs in about 1000 ns ; consequently the input waveform is approxmiately an impulse,

20
which in turn, says that the output waveform will approximate the impulse response of the coaxial cable,

h (k) . Equipped with these physical details of the input and output waveforms, further comments can be

o made about the deconvolution results presented in figures A. 6 and A. 7. Because the output waveform in figure 4.10 was produced by a very short duration input pulse,

variations similar to those in the output waveform should be detectable in the deconvolved waveforms,

d(k) and h(k). Inspection of figure 4.7 shows different oscillations before the main rise of each

waveform. The time scale in figure 4.7 is 50 ns/Div. which is one-quarter of that in figure 4.10. Because figure 4.10 doesn't possess oscillations in the 50 ns time Interval prior to the start of the main waveform rise, those appearing in figure 4.7 are due to noise or the filtering process.

For example, refer to figure 4.6a and note the noise level on the baseline of the waveform; then

look at figure 4.7a which is an expanded time scale version of figure 4.6a. Consequently, it is clear

that the oscillations before the main rise are due to noise.

Next, look at figure 4.7b and note that the regularization has smoothed out the variations due to noise, but that the smoothing very close to the initial rise is not as great as that elsewhere.

Stronger regularization increases the magnitude of the oscillations as evidenced in figure 4.7c. This

is due to the filtering by the regularization process which introduces the Gibbs phenomena due to the truncation of H(n) when forming D(n) [15]. Such effects can be made negligible by avoiding abrupt or rectangular filtering of the H(n) spectrum. Smooth filtering is used in the one-parameter deconvolution

method discussed in Section 4.6.

In Section 4.1 it was pointed out that the FFT and IFFT computations temselves introduce noise or errors due to finite register length. Such errors also cause the imaginary part of the IFFT to be nonzero. For example, say, the real time domain sequence f(k) is given and is then transformed by the FFT to F(n); next, f(k) is recovered from F(n) via the IFFT. However, upon inspection of the data for IFFT {F(n)} it is found that the result is complex with a very small imaginary part while the real part

is very close to the values of the original sequence f (k) . In practice, such non-zero imaginary parts

are always present, and are relatively very small, figure 4.8. In Section 4.5, the nature and origin of the errors affecting d(k) will be discussed.

4.4.3 The Estimation of the Insertion Impulse

Response of a Broadband Antenna

As a second example for the application of the two parameter method, the method was applied to ex- perimental data to obtain the insertion impulse response of a broadband antenna. The data was obtained

using the NBS Time Domain Antenna Range in conjunction with the NBS APMS. Schematically, the experimental arrangement is shown in figure 4.24. Figures 4.25 and 4.26 show the input, x(k), and the output, y(k), and the magnitudes of their

respective DFT's, jx(n)| and |Y(n)|. The magnitude,

|H(n)| = |Y(n)| ^ |x(n)|

is shown in figure 4.27. Deconvolution results are shown in figures 4.28a and 4.28b for two different

degrees of filtering; note that the smoother result and most error occurs for the greater degree of

filtering.

21
4.5 On the Nature of the Error, e(k)

With the results presented in Section 4.4 in mind, it is now appropriate to provide a model for the error e(k) encountered in frequency domain deconvolution. The model is not restricted to any specific

frequency domain regularization filter; consequently, it is valid for both the two-parameter method

(Sections 4.3, 4.4) and the one parameter method (Section 4.6). The discussion will be divided into

three parts:

1. Error in the absence of both regularization (filtering) and computation

errors (noise)

2. Error in the presence of regularization and absence of computation errors.

3. Error in the presence of both regularization and computation errors, being contained in subsections 4.5.1 through 4.5.3, respectively.

4.5.1 Error in the Absence of the Filtering Operation and Computation Noise

The error sequence, e(k), as given by (3.4-1) includes the effects of filtering and computation

errors in the d(k) term,

e(k) = y(k) - x(k) * d(k) (4.5.1-1)

In the absence of filtering and computation noise

d(k) = h(k) (4.5.1-2)

then

e(k) = (4.5.1-3)

because

y(k) = x(k) * h(k) (4.5.1-4)

In words, the error sequence, e(k), will be zero only when d(k) equals h(k); any other choice for d(k)

will give a non-zero result.

4.5.2 Error in the Presence of the Filtering Operation

and in the Absence of Computation Errors

Now consider the nature of h(k) for the signals x(k) and y(k) which are both noisy signals, i.e.,

x(k) = X (k) + n (k) (4.5.2-1)

o X

y(k) = y^(k) + n (k) (4.5.2-2)

where n (k) and n (k) are the noise components associated with the true signals x (k) and y (k)

,

X y o o respectively. The system impulse response is

h(k) = [y^(k) + n (k)](l/*)[x^(k) + n^(k)] (4.5.2-3)

whose DFT is

Y (n) + N (n)

O X

22
The transfer function H(n), (4.5.2-4), obtained from the deconvolution y (k) (l/*)x(k) ) can be

expressed as

1 + N (n)/Y (n)

""^^^ - "o^-) 1 ^ N\n)/X°(n)

- (A.5.2-5)

X o where

Y^(n)

^o^"^ = T~M (4.5.2-6)

Define, °

1 + N (n)/Y (n) ^^"^ =

1 + S^(n)/X^(n) (4.5.2-7)

H (n) is the DFT of the stable true impulse response defined by

h^(k) = y^(k)(l/*)x^(k) (4.5.2-8)

Consequently, H(n) can be viewed as the product of the true transfer function H (n) and some noisy

transfer function, M(n) , which may or may not be stable.

H(n) = H (n) M(n) (4.5.2-9)

o

If H(n) is filtered by R(n) to deemphasize the effect of the noise, M(n) , then D(n) results

D(n) = R(n) [H (n) M(n)] (4.5.2-10)

o

Assuming that R(n) operates only on M(n) , i.e., write

D(n) = H (n) [R(n):M(n)] (4.5.2-11)

o

in the limit

lim D(n) = H (n) (4.5.2-12)

R:M^ 1

°

However, it is physically impossible for the filter R(n) to operate solely on M(n). In practice R(n)

also operates on H (n) with the design of R(n) being such that it affects M(n) more than H (n) and makes d(k) at least stable and possibly smooth.

Consequently, filtering attempts to make h(k) stable by remvoing m(k) from the impulse response,

h(k) = h (k) * m(k) (4.5.2-13)

o

in such a way that

d(k) = h (k) * [r(k) * m(k)] (4.5.2-14)

o

that is, d(k) becomes approximately

d(k) = h (k) (4.5.2-15)

o

23
where the last equation is a judgment by the individual performing the deconvolution. The resultant d(k)

does not cause the error, e(k), to go to zero; only h(k) can do so. Therefore, it can be concluded that a paradoxial situation exists: h(k) makes the error, e(k),

vanish; but h(k) is not the desired result, h (k) , which the individual performing the deconvolution

wishes to achieve. The paradox occurs because y(k) is used as the reference waveform or sequence

in the error expression, (3.4-1). As pointed out in Section 3.4, y (k) is not known; consequently,

the known quantity y(k) is used. Since h(k) and H(n) are given by (4.5.2-3) and (4.5.2-4), respectively, the error, e(k), can be

expressed as

e(k) = y(k) - x(k)*r(k) * [y (k) (l/*)x(k)

]

(4.5.2-16)

= y(k) * [6(k) - r(k)] (4.5.2-17)

where 6(k) is a unit impulse. In the discrete frequency domain, the error is

E(n) = [Y^(n) + Ny(n)] - R(n)

rY (n) + N (n)-] X (n) + N'(n) f^o(-> ^ \^-^^

o X J (4.5.2-18)

which corresponds to (4.5.2-16) while

E(n) = [1 - R(n)] [Y (n) + N (n)

]

(4.5.2-19)

corresponds to (4.5.2-17). The square of the error magnitude is

|E(n) 1^ = il - R(n) 1^ [y (n) + N (n)

1^ (4.5.2-20)

and the error energy over the N points comprising the time window or one period of e(k) is, by the Parseval relation [22, 23],

N-1 N-1 E e^(k) = ^ Y. |E(n)r (4.5.2-21)

where

N k=0 n=0

N-1 N-1 y |E(n)! = y !l - R(n) r |y (n) + N (n)

(4.5.2-22) n=0 n=0 ^

Consequently, it is seen that the error energy per period depends upon the filter characteristics R(n)

Y (n) , and N (n) . The \ o y of (3.4-4), is given byY (n) , and N (n) . The variance or the square of the standard deviation of the error, i.e., the square

o y

, T N-1 _

O = ^ Y, [e(k) - e

]^ (4.5.2-23) k=0

24
? 1^-1? -7

a'' = i y; e (k) - e

" (4.5.2-24)

2 1 ^"1 , ,2-2 a^ = ^ ^ |E(n)r - e (4.5.2-25)

N n-0

where e is the mean value of the error (3.4-3). The effects of filtering on the error can be graphically illustrated by considering the magnitude

of the error, E(n). From (4.5.2-20) the magnitude is

E(n) I = ll - R(n) I |y (n) + N (n)

I

(4.5.2-26)

Figure 4.29 shows the error for two choices for R(n) where R(n) is real; the absolute error, |E(n)I,

approaches zero as R(n) approaches unity.

4.5.3 Error in the Presence of the Filtering Operation and Computation Errors

In the presence of filtering and computation errors the estimated impulse response, d(k), in (3.4-1)

is replaced by d.(k). The error is then given by

e(k) = y(k) - x(k) * d2(k) (4.5.3-1)

which can be explicitly written as

e(k) = y(k) * [6(k) - r(k) - m^(k)] - x(k) * [r(k) - n^^(k) + n^ (k) + n^^{]!.)] + n^(k)

(4.5.3-2)

and which will be derived below. e(k) reduces to (4.5.2-13) when the computation errors vanish, that is,

when

m (k) = 6(k), unit impulse (4.5.3-3)

nj^^(k) = (4.5.3-4)

n^^Ck) = (4.5.3-5)

n (k) =0 (4.5.3-6)

e

The derivation of e(k), (4.5.3-2), and the definitions of (4.5.3-3) through (4.5.3-6) are based upon the assumption that the round-off error in calculating a specific quantity, can be represented

as an additive term to the desired quantity [16]. For example, if A(n) is being calculated from a(k)

by the FFT, then the result is

A^(n) = A(n) + N^(n) (4.5.3-7)

where A(n) is the desired result, and N'(n) is the computation error or noise. Similarly, the IFFT

25
of A, (n) would also introduce an additive error, n . (k)

1 al

a^(k) = a(k) + n (k) + n ^ (k)

1 a al (4.5.3-8)

Thus, it is assumed that the cummulative error depends upon the number of computations, and that, a(k)

is not exactly recoverable from the numerical calculations.

To derive e„(k), (4.5.3-2), the starting point is x(k) and y(k) and the procedure follows below.

FFT [x(k)] = X (n) = X(n) + N' (n) (4.5.3-9)

FFT [y(k)] = Y^(n) = Y(n) + N^(n) (4.5.3-10)

Here, N' (n) and N' (n) are computation errors (noise): note that both X(n) and Y(n) also contain measureX y ment noise corresponding to the noise contained in x(k) and y(k), (4.5.1-5) and (4.5.1-6), respectively.

Next compute the quotient of (4.5.3-10) divided by (4.5.3-9); this result is defined as H (n)

,

H^(n)

Y^(n)

X^(n) + KAn) hi' (4.5.3-11)

where N, ^ (n) is the additive error due to the calculation of H. (n)

hi 1 Filtering H (n) yields D (n)

D^(n) = R(n) Y^ (n) + N^i (n) X (n) hi + N (n) (4.5.3-12)

where N (n) is the additive error due to the calculation of D. (n) . In terms of the theoretical transforms Y(n) and X(n), D (n) is given by

D^(n) = R( ,rY(n) j 1 ^ N;(n)/Y(n)

) "^[^X(n)

i

1 + N^(n)/X(n) ^' ^ ^ (n) + N^(n) (4.5.3-13)

^("^[IS)- >^i(->^\i(")] ^\(-) (4.5.3-14)

where M (n) is a multiplicative error term defined by the { ; term in (4.5.3-13).

The IFFT of D (n) is given by

d2(k) = r(k) "|[y(k)(l/*) x(k)]* m^(k) + ti^-^'^''^ }'^ "r ^''^ "^ ''d2 ^^^ (4.5.3-15)

where n (k) is a complex additive error due to the IFFT calculation of d (k) . The error e(k) is given

by

and

e(k) = y(k) - x(k) * d2(k) + n^(k)

d_(k) = d,(k) + n (k)(l/*)x(k)

5 I e

(4.5.3-16)

(4.5.3-17)

then

e(k) = y(k) - x(k) * d2(k) (4.5.3-18)

26
e(k) = y(k) * [6(k) - r(k) * m^ (k) ] - x(k) * [r(k) * n, . (k) + n^(k) + n,_(k)] + n (k)

.

X n± r dz e (4.5.3-19)

where n (k) is the additive error due to the calculation of e (k) . e When the filtering operation is deleted

r(k) = 6(k) (4.5.3-20)

n (k) = (4.5.3-21)

d_(k) will change to d,(k) and the error becomes

e(k) = y(k) - x(k) * d^(k) (4.5.3-22)

Explicitly, from (4.5.3-19) the error becomes

e(k) = y(k) * [6(k) - m^(k)] - x(k) * l%^(^) + n^2(^)l + n^(^)- (4.5.3-23)

This is the error obtained when no filtering is employed and computation noise is present. When

computation noise is absent, (4.5.3-3) through (4.5.3-6) prevail and e(k) reduces to zero. As pointed out in connection with (4.5.3-15), d„(k) is complex due to the complex additive error

term, n,.(k), which results from the IFFT calculation. Consequently, the estimated impulse response,

dz (4.5.3-15)

d2(k) = r(k) * [y(k)(l/*)x(k) * m^(k) + nj^^(k)] + n^(k) + n^^*^'^)

or d2(k) = r(k) * [h(k) * m^(k) + ti^3^(k)] + n^(k) + n^^^^^) (4.5.3-24)

is complex. Also, the estimated impulse response, without filtering, d,(k), is complex. From

(4.5.3-24) with (4.5.3-20) and (4.5.3-21) d (k) is given by

d^(k) = [h(k) * m^(k) + n^^W] + n^(k) + n^^ (k) (4.5.3-25)

27
4.6 The One Parameter Method

4.6.1 The Filter Structure

This filter is an optimum filter in the sense that it minimizes a performance measure with respect

to the weighted sum of the error power and a smoothness power [18,19]. It turns out that there is only

one adjustable parameter in this method, and the parameter adjusts the weighting between the error and smoothness. The derivation of the filter is given below. Given the two signals x(k) and y(k)

x(k) = X (k) + n (k) (4.6.1-1)

y(k) = y^(k) + n (k) (4.6.1-2)

computation and filtering lead to the estimated impulse response, d„(k), (4.5.3-16),

d2(k) = r(k) * [y(k)(l/*) x(k) * m^(k) +n^^(k)] + n^ (k) +n^2(^) (4.6.1-3)

and the error function.

e(k) = y(k) - x(k) * d2(k) (4.6.1-4)

In addition to the error constraint, e (k)

define an auxiliary constraint, s(k),

s(k) = c(k) * d2(k) (4.6.1-5)

where c(k) is the discrete function that defines the constraint. In the present discrete case c(k) is the second-difference operator, which corresponds to the second derivative operator for continuous functions,

2 2 d /dt , thus s(k) depends upon the rate-of-change of the first derivative of d. (k) , i.e., its smoothness.

N-1

E = XI [e(k)] , (4.6.1-6)

k=0

_ N-1

S = E [s(k)]'' (4.6.1-7)

k=0

Now, to obtain an optimum solution, adjust d (k) until the performance measure, P,

P = E + Y s" (4.6.1-8)

is minimized. Y is a weighting function which selects the degree of smoothing. Consequently, the minimization condition is approximately expressed by

3d2(k)

E + Y S

I (4.6.1-9)

28
assuming that N is a very large number so that d-(k) approximates a continuous function rather than a discrete one.

In the DF domain d (k) , e (k) , and s(k) become

°2(") = ^(">[lSj

M^(n) + N^^(n)

E(n) = Y(n) - X(n)D2(n)

S(n) = C(n) D2(n)

+ N (n) + N^^(n)

r d^ (4.6.1-10)

(4.6.1-11)

(4.6.1-12)

Parseval's relation [22,23]

E a (k) k=l

1 N-1

n=0

(4.6.1-13)

jives for E and S,

N-1

n=0

E(n) E (n) (4.5.1-14)

N-1 V S(n) S"(n)

n=0

(4.5.1-15)

where E (n) denotes the complex conjugate of E(n), etc. The sum E + Y S is given by

N-1

^ + "^ ^ " N ^ ^^^^^ E (n) + Y

S(n) S (n) ] n-0

(4.5.1-16)

where

E(n) E (n) = [Y(n) - X(n) D2(n)] [Y(n) - X(n) D2 (n)(4.6.1-17)

S(n) S (n) = C(n) D2(n) C (n) D2(n)

|C(n) D2(n) (4.6.1-18)

Consider E(n) E (n)

E(n) E (n) = [Y(n) - X(n) D2(n)] [Y(n) - X(n) D2(n)]

[Y(n) - X(n) X)^M] [Y"(n) - {X(n) D2(n)}"]

= Y(n) Y"(n) - Y(n) {X(n)D^(n)}" - Y"(n) {X(n)D (n)}

+ {X(n) D2(n)} {X(n) D2(n)}

29
E(n)E*(n) = |Y(n)|^ - Y(n) {X(n) D2(n)}* - Y*n) {X(n) D^Cn)} + |x(n) D(n)|^

(4.6.1-19) Hence

E + Y S = I 2 I |Y(n) 1^ + |x(n) D^Cn) |^ + Yic(n) D^Cn) |^ - Y^di) {X(n) D^di)}

n-0 '

- Y(n) {X(n) D^Cn)}*! (4.6.1-20)

Consider X(n) D (n) . From, (4.6.1-10),

X(n) D^dn) = R(n) [Y(n) M^(n) + X(n) Nj^^(n)] + X(n) [N^(n) +N^2(n)] (4.6.1-21)

First, Im. (n) I is of the form 1 + J(n), where J(n) is complex and |j(n)|<<l. Second, |n, ,(n)l, In (n) I

1 hi r and |N,„(n)| are the magnitudes of computation errors; consequently, they too will be much less than

unity. When X(n) D„(n) is written in the form

X(n) D2(n) = R(n) Y(n) A(n) (4.6.1-22)

A(n) will be complex and of the form

A(n) = 1 + K(n) (4.6.1-23)

where

|K(n)| « 1 . (4.6.1-24)

Consequently, A(n) will be very close to the real number unity.

Some complex relations:

XD = RYA (4.6.1-25)

= IrIIyIIaI * +4 + d) (4.6.1-26)

. r y a

(XD)* = IrI IyI IaI /-((t> +4) + <t> ) (4.6.1-27)

^ ^ I I I I I

1 / ^r y __a

Y = |y| /4) (4.6.1-28)

Y YI/-4 (4.6.1-29) ^1

y""'(XD) = IrMyI lAJ/cf) +4) (4.6.1-30)

Y (XD)* = 1r1 IrMyI^IaI |y1^|a| /-((}) J-{<p^ + c}) ) (4.6.1-31)

Y*(XD) + Y(XD)""' = |r1|y|^|a| {2 Cos (<[) - <f )

(4.6.1-32)

r a

= 2R|y|^!a| Cos (j) (4.6.1-33)

where

Putting (4.6.1-33) into (4.6.1-20) yields

(4.6.1-34)

30
N-1 E+YS = I X! ||Y(n)|^ + |R(n) Y(n) A(n)I^ + Y|C(n) D^CnXl^

n=0

-2R(n) Y(n) | A(n) Cos (})( »>} (4.5.1-35)

i

I;'

|Y(n)|2-. |R(n)Y(n)A(n)i2+Y|^(^i^y^(^^^-A(nl

n=o

- 2R(n) Y(n) A(n) Cos 0( »') (4.6.1-36)

i E I^WlM l^R^n) |A(n)|^+ ^ | C(n) | ^R^(n) |A(n) | ^

" --'^ ' n=0 |x(n)r

2 R(n) |A(n) ] Cos <t>in}\ (4.6.1-37)

In principle, the performance measure (4.6.1-8)

P = E + Y S"

can be expressed in terms of discrete time domain quantities or in terms of discrete frequency domain

quantities by virtue of the Parseval relation (4.5.2-21). However, in practice, given the time domain

quantities x(k) and y(k) , subsequent numerical computations introduce errors into the frequency domain

quantities. Consequently, in the frequency domain the result is ?',

P' = E' + Y

^'

which is not equal to P. To minimize P, d (k) is adjusted. Varying d (k) corresponds to adjusting the frequency domain filter R(n) . Consequently, P and P' are minimized through the relations,

9 P ^ d2(k)

3 P'

a R(n)

=

respectively. Because the computation errors are very small, a first order approximation is to assume

that P is ec R(n), i.e..that P is equal to P'. With that assumption being made a change in d.(k) corresponds to a change in

and

3d2(k) [^ + ^0= 3fao["^^']

-, ^, ,

i E + Y S I = .i I, .] E + Y sl + 1^-7- E + Y S

9 R(n) {_ 'J 3 1 R(n) 1 L J 3 <t>_ L J

(4.6.1-38)

(4.5.1-39)

Consequently,

31
C^[e + Ys] = i E |Y(n)|2 (2iR(n)i!A(n)|2 +

2y | C(n)A(n) | ^

J RCn) n=0 |X(n)

- 2|A(n) |Cos[(})^(n) + (f)^(n)] + 2 | A(n)R(n) | Sin[(t)^(n) + (J)^ (n) ] i (A. 6. 1-40)

To minimize [E + Y S] with respect to R(n) , lR(n)| and (f) (n) must be such that (4.6.1-40) is zero.

I 1

2

Since |Y(n)] is not generally zero, the { } -term in (4.6.1-40) must be made equal to zero to force

(4.6.1-40) to zero. Choosing the filter phase angle to be the negative of the computation error phase

angle <t>^M .

()) (n) = -<p (n) (4.6.1-41)

eliminates the phase angle dependence, and gives for |r (n) | ,

|R^(n)| = ^^^ T (4.6.1-42)

lA(n) I |x(n) I + Y|A(n) | |c(n) \

Consequently, the optimum filter |r (n) | /(}) (n) is given by (4.6.1-41) and (4.6.1-42) in the presence

of computation noise.

If smoothing is not employed, 7 is set equal to zero; and the filter becomes

|R (n) I = u/. I (4.6.1-43)

' o 'y=o |A(n)

I

[*ro^")lY=0 ^ *a^"^ (4.6.1-44)

which says that the filter magnitude and phase are such that they just cancel out the computation errors

so as to force the error function e(k) to zero. Furthermore, when computation errors are not present

A(n) is unity, i.e.,

|A(n)| = 1 (4.6.1-45)

cf) =0 (4.6.1-46)

which, in turn, means that (4.6.1-43) and (4.6.1-44) reduce to

[R(n)] = 1 (4.6.1-47) Y=0 A(n)=l

[<})„(n)] = (4.6.1-48) ^ Y=0

A(n)=l

which is the case of no filtering, i.e., R(n) = 1. As R(n) approaches unity, the filtering l-R(n)

approaches zero, figure 4.29.

32
In the presence of smoothing, Y ^ 0, and in the absence of the computation errors, A(n) = 1, the optimum filter is real and is given by

I

1

R(n) = ^^^^ (A. 6. 1-49)

|x(n)r + 7|C(n)r

Before discussing the form of C(n), some summary comments are in order to clarify the meaning of the optimal result in the presence of computation errors, (4.6.1-41) and (4.6.1-42). It must be kept in mind that in any given problem A(n) is not known, i.e., |A(n)| and (\> (n) are unknown; consequently, the optimal filter cannot be calculated. Equations (4.6.1-41) and (4.6.1-42) are only of theoretical

interest in that they show that computation errors require that the optimum filter, R (n) , be different

from that in the absence of computation errors, R(n) . Specifically, R (n) equals A(n) R(n) . In practice,

since A(n) is not known, it is convenient to let R(n) be real and given by (4.6.1-49), which is the case

in the absence of computation errors. Again, because A(n) is very close to the real value of unity, the use of (4.6.1-49) is the best that can be done.

c(k) is defined as the second difference sequence operator and corresponds to the second derivative

operator in the continuous case. The first backward-difference sequence, V, is defined by the operation

{x(kT) - x[k-l)T]}/T where the interval, T, between the k-th and (k-l)-th values is explicitly shown. Implicitly, the operation is denoted as x(k) - x(k-l). V * x(k) is given by

V * x(k) = [1,-1,0,..., N-1] * x(k) (4.6.1-50)

2 The second backward-difference sequence, V , is defined by the implicit operation x(k)-2x(k-l)+x(k-2)

2 V "' x(k) is given by

V^ * x(k) = [1,-2,1,0,..., N-1] (4.6.1-51)

2 Note that there are N-3 and N-4 zeros in the V and V sequences, respectively. The DFT of

V^ = [1, -2, 1, 0, 0, 0, N-1.] (4.6.1-52)

is given by the complex sequence

A^ = ]^ V^(k)e ^"^'^'^, n = 0,1,2,0,0,0,...N-1 (4.6.1-53)

k=0

1 - 2e"^^ '^^ + e"J^ "^"^

; n = 0, 1, 2,0,0, 0, .. .N-1 (4.6.1-54)

which also contains N-4 zeros, and where

n T = -i^ (4.6.1-55) N

2 Since c(k) is defined as V (k) , then C(n) is given by

. 2 TTn . 4 TTn C(n) = l-2e + e , n= 0, 1 ,

,0,0,0, . . .N-1 (4.5.1-56)

33
and

C(n) = C(n)C(n)

. 2 TTn . 4 Tin

, "^ N _^ "-' N l-2e + e

. 2 TTn . 4 TTn

, . ^ N , ^ N l-2e + e

n = 1,1,2,0,0,0,.. .N-1 (4.5.1-57)

|c(n)|^ = 6 -8 Cos ^-^ + 2 Cos ^Y^ ; n = 0,1, 2,0,0,0, .. .N-1

which can be reduced to [21]

. 4 C(n) = Sin (2TTn/N); n = 0,1, 2 , 0,0,0, . . . ,N-1

(4.6.1-58)

(4.6.1-59)

When N becomes very large, the argument, 2TTn/N, approaches a continuous function of f; also, the argument becomes very small, i.e..

large

N

L N J

2 TTf, <f< 1. (4.6.1-60)

Since

Sin oj = oj- —^ +

3 5 7

0) . OJ w

' 7 ' 1\ + • • •

; oj < CO (4.6.1-61)

then for

^^^^^ [C(n)|^ = |c(jw)|^ = 0)^ (4.6.1-62)

which corresponds to the square of the magnitude of the continuous function second derivative operator.

2 2 D^ = (j oj)^ (4.6.1-63)

dt

D = w (4.6.1-64)

34
4.6.2 Conversion of a Step-Like Waveform into a Duration Limited One

Given a waveform or sequence, f(k), comprised of N values, successive application of the DFT and

IDFT yield an F(j(jj) and f (k) which are periodic of period, N. Two examples are shown in figure 4-30,

where the given sequence of N points are at the left-hand side of the figure. The resultant IDFT's are periodic in which one cycle is identical to the original f(k). For the sequence f (k) , the spectrum F (n) corresponds to the intended-time sequence f (k) . In

the case of f, (k) the spectrum G (n) does not correspond to the intended-time sequence. Consider the

b b following physical situation, the step-like waveform f, (k) is N points of a sequence that begins at zero and continues to n = >», remaining at the constant value B for k J^ k. where k £ N-1. Because the DFT will use only N points, the f, (k) is abruptly truncated at k = N-1. The resultant spectrum G, (n) is

b b not the intended-spectrum, F (n) . The IDFT shows g^ (k) , a cycle of which is not the intended f, (k)

b b b i.e., the constant step-like waveform. The abrupt truncation f. (k) introduces the spectral components

b of an abrupt transition (step function). A step-like sequence, f(k), which has achieved a constant value, B, for k. £ N-1, can be converted to a duration-limited one, of length M = 2N, figure 4.31. The conversion is accomplished by extending and delaying f (k) by N units to obtain f(k-N), figure 4.31a and b. f(k) and f(k-N) both consist of M = 2N points. f(k-N) is subtracted from f (k) to yield f„(k).

f (k) = f(k)-f(k-N), k = 0,1,2, •••M -1 (4.6.2-1)

The DFT of f„(k) is given by

F2(m) = F(m) - F(m)e J^^TM^, m = 0,1,2, •••, M-1 (4.6.2-2)

2tt where F(m) is the DFT of f(k) and Q = — . F (m) is the product

where

Because

F^Cm) = (1 - e i^'^^^) F(m) m = 0,1, 2 , • • • ,M-1 (4.6.2-3)

F (m) = F (m) F(m) (4.6.2-4)

F^(m) = 1 - e-^"^"^ = 1 - e-^'"^ m = 0,1,2, •••, M-1 (4.6.2-5)

|f (m)| = 2 |Sin mTT/2| (4.6.2-6)

L(m) = Tan-1 TT^^ (4.6.2-7)

1 1-Cos mTT

It^/nI 0,m=0,2,4--* //coQ\ I F (m) = (4.6. 2-0)

1 2, m = 1,3,5-

|)^(m) « 0, m = 0,1,2, ••• (4.6.2-9)

35
the DFT's of F(n) and F (m) are simply related. Figure 4-32 shows F(n) and F (m) where m has been used

in the argument of F . Note that [ F (m) j is twice the value of |F(n)| for m = 2n-l. Explicitly,

F(n) and F (m) are F(ni;) and F [m(n/2) ] . The discrete frequency nfi corresponding to n = N-1 is equal to the discrete frequency mfi/2 with m = M-2 because M = 2N. Thus, the spectrum of F (m) has (a) the same number of non-zero values, N, as that for F(n), but are graphically spaced twice as far apart, and (b) the'magnitude of JF.(m)[ is twice that of |F(nn)| at the non-zero values.

In deconvolution studies the effect of the time domain conversion of step-like waveforms for x(k) and y(k) yields an estimated impulse response due to a unit impulse applied at k=0, followed by a negative unit impulse applied at k=N-l, figure A. 33.

36
A. 6. 3 Application of the One Parameter Method

4.6.3.1 Procedural Sunimary

In the discussion presented here, the notation will be that employed in Section 4.5; consequently,

d(k;) will be more precisely specified by d„(k); similarly D^ (n) will replace D(n), and so on for other

quantities as the occasion arises. The reason for doing so, is to emphasize that the quantities contain

errors. In the sections before 4.5, the errors were not delineated.

The procedure starting with the time domain sequences x(k) and y(k) and ending with the deconvolved

time domain impulse sequence d (k) progresses as follows:

1. Offset correction for x(k) and yCk). Eliminate the offset, i.e., insure

that x(o) = x(N-l) = and y(o) = y(N-l) = 0. The offset error is systematically introduced by the data acquisition device. If x(k) and y(k) are step-like waveforms so that x(N-l) and y(N-l) are not zero, they must

be transformed to impulsive waveforms. See subsection 4.6.2.

2. DF transform x(k) and y(k) to obtain X,11(n) and Y (n) ; display [x (n)l and

1 do |Y^(n)|^g, (4.5.3-9), (4.5.3-10).

3. Calculate the complex ratio H Ill fn) = Y^(n)/X^(n), (4.5.3-11); display |h fn)i ,^. 1 dB

4. Select the regularization filter coefficient y; then filter H (n) to yield

D^(n), (4.5.3-14) display |D^(n)|^g.

5. IDF transform D (n) to obtain the time domain deconvolved result,

d2(k). (4.5.3-14); display d (k)

.

5. Calculate the error criteria (4.5.13-19); display e(k) and the values for e , _ max.

e . , e, and 0. mm.

7. Judge the quality of the deconvolved impulse response estimate, d (k) , in terms of the error criteria; if unsatisfactory, return to step 4 and change

the filter parameter.

8. If the quality is satisfactory, print the d (k) graphic display, related error

criteria, and the numerical values.

Three examples of the one-parameter deconvolution method will now be presented; they are:

1. The estimation of the insertion impulse response of 300 meters of RG 58 C/U coaxial

cable; this example uses the same insertion waveform data as was used in the two- parameter method example. Section 4.4.1.

2. The insertion impulse response of a filter for a step-like, 100 ps transition

duration, reference waveform generator.

The first example allows the reader to compare the one and two parameter methods; the second

example, in addition to illustrating the one parameter method, illustrates the procedure for using

step-like waveform data rather than duration limited waveforms.

37
4.6.3.2 Estimation of the Insertion Impulse Response

of 300 Meters of RG 58 C/U Coaxial Cable

The insertion waveforms x(k) and y(k) are the same ones used in the two parameter method experiments

in Section 4.4.2, figures 4.9 and 4.10, respectively. Three degrees of regularization were used, Y

= 100,

500, and 10,000. The estimated insertion impulse responses and their corresponding error functions, e(k),

and filtered (regularized) spectrum magnitudes are shown in figures 4.34 through 4.36. The regularization

filter was that of (4.6.1-49),

„, , |X(n) 1^ R(n) = ' %' ' r lx(n)r + Y|c(n)r

which is a real filter.

When Y

= 100, the regularization is moderate. Compare figure 4.34 with figures 4.7b and 4.12 of the two parameter method. The results are comparable. Increasing Y to 500 provides strong regularization. Compare figure 4.35 with figures 4.7c and 4.13 of the two parameter method; note that the Gibbs

phenomena has appeared in the results of the two parameter method, figure 4.7c, while the one parameter

method produces no Gibbs oscillations. Increasing Y to 10,000 provides still stronger regularization

with no evidence of Gibbs oscillations (compare figures 4.36 to 4.7c and 4.13). Furthermore, figure

4.36 very clearly shows the small variations (reflections) discussed in the sixth paragraph of section

4.4.2. Thus, it is seen that the one parameter method uses a smooth filtering that does not introduce

the Gibbs effect (as was pointed out earlier in section 4.4.2).

38
4.6.3.3 The Insertion Impulse Response of a Filter for a Step-like,

100 ps Transition Duration, Reference Waveform Generator

In this example, the insertion impulse response of a liquid filled coaxial transmission line is estimated. The coaxial line is about 20 cm long and is filled with a dilute binary solution in which

heptance is the solvent and heptanone is the solute. Because heptanone Is a polar molecule (heptane

is non-polar), the solution comprises a dilute polar dielectric and as such, exhibits a dielectric

relaxation law. Consequently, the frequency dependence of the dielectric constant, e(s) causes the transmission line to behave as a low pass filter. The filter is used to shape a 20 ps transition

duration (10%-90%) step-like waveform into a 100 ps transition duration waveform [25].

The insertion step-like waveforms x(k) and y(k) are converted to duration limited waveforms, X-. (k) and y^ (k) , and the DFT's X (n) and Y (n) are obtained, respectively, figure 4.37. The effects of filtering on the DFT's are shown in figure 4.38; here, the real filter, R(n), of (4.6.1-49),

Rfn^ - |X(n)|^

R(n) - y' r

ix(n)|^ + 7|c(n)r

2 4 6 — —

has been used for 7 = 0, 10 , 10 , and 10 . R(n) operates on H (n) to yield D (n)

D^(n) = R(n) H^(n)

and the IDFT of D (n) is d„(k). The estimated impulse response for very strong regularization, Y

= 10

is shown in figure 4.39. Note that d (k)

d2(k) = y(k)(l/*) x(k)

is the response due to a pair of impulses as shown in figure 4.33, and d (k) is the desired estimated

impulse response, the first pulse of the pair in d„(k).

This example completes the discussion of frequency domain deconvolution. In the next chapter,

some time domain deconvolution methods will be presented.

39
K^(k) y„(k)

Figure 4.1 X (k) and y (k) are the exact sequences which are uniquely related through the system impulse

response sequence h (k)

x^(k)

Measurement

System

t

x(k)

h_,(k)

1

1

n^(k)

y„(k)

Measurement

System

t

y(k)

n^(k)

Figure 4.2 x(k) and y(k) are corrupted versions of x (k) and y (k)

The noise sources n (k) and n (k) are independent of

X y

each other.

40
d^(k)

d2(k)

d3(k)

d^(k)

^* k

|-^- k

n-1

Figure 4.3 The effects of filtering, d (k) , the stable response resulting

from the "proper" amount of filtering, h (k) in the operator's

judgment, d (k) , almost stable, needs more filtering. d^Ck),

very unstable, very little filtering, if any. d (k) , stable,

too much filtering, response is slowed down.

41
dB

H(no)

100 dB

o o Discrete Frequency, nf.;

Figure 4.4 The magnitude, |D(n)|j_,, of the filtered H(n) versus n Q, dB

in dB, 20 log |D(n) | , and the magnitude of the unfiltered

|H(n)1^3.

42
50 fi

e(k) 50 fi MEASUREMENT

SYSTEM

x(k) = X (k) + n (k)

O X

50 fi —AAASYSTEM

UNDER

TEST

MEASUREMENT

e(k)^ SYSTEM

mmmm -y-y w^

50 ^.^ y(k)=y^(k)+ny(k)

Figure 4.5 Insertion response measurement method. The observables are the

waveforms x(k) and y(k) obtained before and after the system under

test is inserted between the 50 9, impedance generator and

measurement system, respectively. H(n) = Y(n)/X(n) = S_ , the

insertion scattering parameter [17].

43
RE«L n«T CF IX OECONUCLUriON RESULT KM. PWT OF TtC OECOMUUTION ICSULT

I

ni« • e.zzx -1 \,

N. niN —e.25BE -2 s JIN

, H, m^11

HORIZOHTM. SCALE 1

1 r

1

DIU Zn.ei NB I\

^

^ s^•*• k^

i^k i_ MP KPf «.». r

DIU • e ZSSE -2

NPK - e.22SE -1

NIN a-eZSK -2

1

HDRIZOHrM. 90MX In [

DIU - 298.01 HB ,1 m ji,

llf L, V^

N

. . .iiiinlil :r.yr^'i,.t

riTl

III I.I I. II l.ii m

•'

I I I "Ml ! ! •;

FUX SPECTRUH

(a)

KfiL PMtT OF T>C OECONUOLUTIOH RESULT

DIU • e 2SK -2 V WW • a.22SE -1 \

niN a-e.TKF -2 \-*^ Li

s iJ M

1 **! yj At In

HORIZOHTM. SCfiZ r 1 i! }} M

DIU - 2ee.ei he \

-l

\ i Hki

1 Ml"Hi •A N. ku^U<J

n

^

"TW^n^ '•W

^^^^H 1 ,1 FULL SPECTKUn

(b)

REAL PART OF THE OECCMnU/TION RESULT

> J —t

T| 1 iJ

riiN ^-e.zseE -2 Nt\

1. 1 ,ii y

i.'T(i,i^n::,i m

1, '^raiPEi 1

HQRIZCHTM. SOLE Vi

1

DIU • aaa.ei te ri

L . ii:VllJilUlklLiMiliL JLiil[i>LliU

iUHIKiiliyifaAkiiiiiikkyiLib

FULL SFGCTRUn FULL SPECnSfl

(C) (d)

Figure 4.6 Deconvolution without any regularization, i.e., full spectrum of H(n) is

used; h(k) is shown for various input and output signal to noise ratios.

(a) SNR-x 75 dB, SNR-y 46 dB. (b) SNR-x 75 dB, SNR-y 36 dB. (c) SNR-x

45 dB, SNR-y 46 dB. (d) SNR-x 40 dB, SNR-y 40 dB. |H(n)L^ is shown in

db

the upper right hand portion of each figure with the vertical scale of

-10 dB/division commencing at 0.0 dB and the horizontal scale of

25 MHz/division commencing at 0.0 MHz.

44

I
A 1 A

1 1 1

\ \

\ \

V

V

\^ \ \ X

H

^"^^ "^

(a) (b)

|i

\

V\

^.

AA A ^-^ r^y

A

\

\

V

\\

"^

h

\

\

\

\\\ X

j-laA

"^

(c)

A

\

\

\

I

\, ^

A

^—\>

Am

(d) (e) (f)

Figure 4.7 Deconvolution with regularization; d(k) is shown, (b-f), for various

input and output signal to noise ratios and degrees of filtering.

SNR-x 75 dB, SNR-y 46 dB : (a) h(k) , no regularization, n ^ = 250 MHz,

A = 0; (b) moderate regularization, n S7 = 100 MHz, A = 5; (c) strong

regularization, n fi = 74 MHz , A = 6. (d) SNR-x 75 dB, SNR-y 36 dB,

strong regularization, n SI = 74 MHz, A = 6. (e) SNR-x 45 dB, SNR-y

46 dB, moderate regularization, n Q = 100 MHz, A = 5. (f) SNR-x

40 dB, SNR-y 40 dB, strong regularization, n U = 74 MHz, A = 6.

Vertical scale 2.5 x 10~-^/Div. , Horizontal scale 50 ns/Div.

45
mm. rmr op iw oBONuajunoN «iut

OIU • • ZSBE -2 mx - e 22S -1

I

fl.5 \,

niN •^.ZSBE -2 \,

s

H,

HDRIZDKTM. 90AL£ 1

1 'v. ^

oiu • an.et w

V\ >^^

OOENT OF T»C U6EFU. aPECTRm • ISe »K

EWOR FIWCTION

- 8 478E -2

• e 243E -1 —«.233E -1

- e ie6E-5

P' a.483E -2

1

. . 1 .III. niN

810*

fr Mimiiiiniii'jfi'iNrfffinr

'' fir '"niiTTi'iii • iin"-'-!- 1 1

1 iMttiiiiluiJui,ii!iii:u>y.

1 T 1 '

ExrE^a cf nc useful vccnu* • in »c

nAUNMlY MRT OF T>E OECaWJOLUTlOH KajLT

EXTEHT CF T>E UKFUL VECT1U1 > 1« »C

Figure 4.8 An example of the imaginary part of the deconvolution result.

If the computation procedures of the DFT (FFT) were exact,

the imaginary part would vanish.

46
INPUT iMUEFORn

oiv; - e.reeE i mx - e.€3eE 2

HQRIZQHTA. SCALE

oiv > 2ee.ee hs

SIOML TO NOISE

7S.1 DB

_^— ^B.—_ ^i_^ _^^_ ^.^^ •_»• i_^ _^_

' __.—^— _^^_ .^_^_ ^^__ .^^^ _^^ ^.^^

rVJOfllUDE OF Tf€ I^fVT MWJEFORM FOURIER TRANSFORM

le O&^IU

HORIZOKTAL SCALE

oiu « 25.ee wz

Figure 4.9 The input waveform x(k) and the magnitude of its

DFT, |x(n)|. SNR-x 75 dB. x(k) is shown with

an abscissa of 50 ns/div. in figure 5.15.

47
OUTPUT WAUEFORM

r--"'

Diu * e.seeE o mK = 0.45eE 1

MiN —e.seeE e

HORIZONTAL SCALE

DIU « ase.ee ns

SIOM. TO NOISE

46.5 OB

tVQ^UUCE CF THE OUTPUT MAUEFORM FOURIER TRANSFORH

10 DS/^IU

fCRIZOKTAL SCALE

DIU » 25.ee MHZ V\ s\

\\1 .1 1 • 1

1

1 it7i:ri"'^p^."n '^v

I'l T"^ 1^ 'f

• fPM V

1

Figure 4.10 The output waveform, y(k), and the magnitude of

its DFT, |Y(n)|. SNR-y 46 dB. y(k) is shown

with an abscissa of 50 ns/div. in figure 5.15.

48
REM. PORT CF T>C OECOKUOLUTION RESLLT

DIU » 0.2S0E -2

nWi » 9.225E -I

WIN —«.25eE -2

HORIZOHTAL SCALE

DIU > 288.91 NS

k

\j> \,>

s lUlU

, H|

If! 1

1 FP

\

V

fM

\

s,Mk N% NH» Hm Mh

FULL SPECTRUh

ESFOR FUNCTION

FULL SPECTRUM

Figure 4.11 SNR-x 75 dB, SNR-y 46 dB. Deconvolution without regularization, lH(n)|^g,

h(k), and e(k). |H(n)| ,„ is in the upper right hand portion of the upper

dB

figure; vertical scale -10 dB/Div. commencing at 0.0 dB, horizontal scale

25 MHz/Div. commencing at 0.0 MHz. No regularization used; this relatively

smooth result for h(k) is possible only because of the large values of

SNR-x and SNR-y (in particular SNR-y) . Note the small value of the error

function e(k) and that its shape is comparable to y(k).

49
REPL P(«ir OF TIC CECONUOLUTION RESULT

oiy e.?rwF -2

rtw « e ??SF -1

V

B-5

\

MIN —«.250E -2 \,^s

i H|

HORIZDKTAL SCALE

1 x

\.

OIU * 268. ei NS \

V\

EXTENT OF THE USEFUL SPECTRUM 188 tfC

ERROR FUNC7ICW

ji jiiiljUjiiiijMiiijfliikyi

1

' 1 1 1

EXTENT OF HE. USEFUL SPECTRUH - 188 mZ

Figure 4.12 SNR-x 75 dB, SNR-y 46 dB. Deconvolution with moderate regularization

n f^ = 100 MHz, A = 5, o |D(n)L„, d(k), and e(k). Note the filtering

dB

D(n) j^ is shown in the upper right hand

dB of H(n) to obtain D(n),

portion of the upper figure, Vertical - 10 dB/Div. commencing at

0.0 dB, horizontal 25 MHz/Div. commencing at 0.0 MHz.

50
RE«. PftRT OF T»C OECONUOLUTION RESULT

Diy O.250E -2

mx « 225E -1

MIN —«.250E -2

HORIZONTAL SOM^

OIU « 288. ei NS -1

I

/9=6

\> \,^N

\s s. •V

\. \

-TrrJ

EXTENT OF "DC USEFUL

n

SPECTRUM - E!WR FUtCTlCW

DIU « 0.227E -1

Tw: « 112E e

PIN —e.llSE

^EPN > e 106E -5

SIOIO> 0.114E -1

74 mz

ff

m. .jii

^^HiIM^UiyyimIHm

m

m

m

EXTENT m

OF T« USEFUL SPECTRUM - 74 WZ

Figure 4.13 SNR-x 75 dB, SNR-y 46 dB. Deconvolution with strong regularization

n ^ = 74 MHz, A = 6. |D(n)|, d(k), and e(k). Note the filtering of

o H(n) to obtain D(n). |D(n)| is shown in the upper right hand

do

portion of the upper figure, vertical -10 dB/Div. commencing at

0.0 dB, horizontal 25 MHz/Div. commencing at 0.0 MHz.

51
QtlTPUT UM.'ER]f91

Diu - e.seeE e mx « 0.45eE 1

«iN —e.seeE e

HORIZQKTAL SCALE

oiu « 29e.ee hs

SIQM. TO NOISE

36.5 oe

\

i V

i >Sum,.m-»Mi -m,

flAWITire OF T>€ OUTPUT WAUEFORM FOURIER TRflNSRJRM

le DB/DIU

[

V

HORIZONTAL SCALE \1^

OIU = 25.ee mz Nll.i. .Hill hi til .III. ihiii ii .1 II II .1

J

mill

1

yUi^.iiiiui.aLiiUiiL.i.1

'iri t " 1 1 i|

r i

Figure 4.14 The modified output waveform, y(k), and the magnitude of its

a SNR-y of 36 dB

DFT, ,Y(n) I , after the addition of 10 dB of noise to yield

52
REAL PPRT OF THE DECONWDLUTIOH RESU.T

DIU - e.258E -2

HAX » e.22SE -1

HIN —«.258E -2

HORIZONTAL SCALE oiu = 2ee.ei hs

FULL SPECTRUM

ERROR FUNCTION

Diy - 3.441E -4

rttX « 9.390E -3

«IN »-O.S10E -4

rcAN « e.ie2E -«

SiaiA= 8.443E -4

I I . ll ll irrriiii

HMMM

J

I

ll hill.

I \

HP|

V

JF>^ m0*^ rW

FULL SPECTRUM

Figure 4.15 SNR-x dB, SNR-y 36 dB. Deconvolutlon without regularlzatlon _ 1h('^)1j3> h(k) and e(k) [H(n)|jg Is In the upper right

hand portion of the upper figure, vertical scale -10 dB/Dlv.

commencing at 0.0 dB, horizontal scale 25 MHz/Div.

commencing at 0.0 MHz.

53
PEW. PPKT OF THE DECONUOLUTIOH RESULT

DIU e.2saE -2

rwx » e.22SE -J

fllN »-«.250E -2

HORIZONTAL SCM^

DIU « 28e.81 NS

*L fl'.6 \

^

1 1 uTli| 'V^

V.^

i ^^

\

\

\ -

\

'-v, ^fySn **¥¥! "r^i"

EXTENT OF THE USEFUL SPECTRUH

ERROR FU^CTION

74 MHZ

DIU e.317E -1

nAX = 8.132E e

MIN *<.164E e

rcfiH « e.ieiE -5

SIQW> e.SlSE -1

.l.t 1 1 I .1.

•'

ll'';[n''^ fill' ^ i f|[';i^ 1

^ • in I

I I ill.i^l.lli ilUlli:ii.ili,Jlki.

I

iq[]i I

III' V(!!

EXTENT OF THE USEFUL SPECTRUM 74 t*C

Figure 4.16 SNR-x 75 dB, SNR-y 36 dB. Deconvolution with strong regularization,

n ^ = 74 MHz, A = 6. o of H(n) to obtain D(n).

portion of the upper figure. Vertical -10 dB/Div. commencing at 0.0 dB,

horizontal 25 MHz/Div. commencing at 0.0 MHz.

54

iD(n)|,„, d(k) and e(k). Note the filtering

' dB

|D(n)| is shown in the upper right hand

dB
INPUT UrtUEFORn

ciu = e.TecE I

nftX = 9.63aE 2

HORIZONTAL SCALE

DIU * 268.08 NS

SIGNAL TO NOISE

45.1 oe

nWINITUOE OF THE IMVT WAUEFORM FOURIER TRANSFORM

18 08/OIU

HORIZONTAL SCALE

DIU » 25.88 mZ

.,-. ^*-^*N "-^^ Will 1 **\J ^• 1. II .1

1

1

1

hmrmr^' m 1

r !•

-

Figure 4.17 The modified input waveform, x(k) and the magnitude of its DFT,

|X(n)L , after the addition of 30 dB of noise to yield a dB

SNR-x of 45 dB.

55
»«ftt, P««T OF T>€ DECONUOLin'IOH REStJLT

DIU - 0.250E -2

fWX « e.225E -1

niN —e.250e -2

HORIZOHTAL SCALE

DIU * 286. ei HS

FU-L SPECTRUM

ERROR RUCTION

r

5 \ -^B N iii

-1

HfyMffn t fWTr^

t

t 1 M "^^»^h.^N,^

FULL SPECTRUM

Figure 4.18 SNR-x, 45 dB, SNR-y, 46 dB. Deconvolution without regularization.

H(n) l^g, h(k), and e(k) |H(n) ! is in the upper right hand

dB

portion of the upper figure, vertical scale -10 dB/Div.

commencing at 0.0 dB, horizontal scale 25 MHz/Div. commencing

at 0.0 MHz. ..
f^A^ Pl«T OF THE OECONUQLUTION RESULT

DIU « e.250E -2

f»X = 9.225E -1

/niN »-e.250E -2

HQRIZONTM. SCALE

OIU - 288.81 NS -\

\

\

f)::5 \V

s

H.

1 ^

-s^

"^

^ f ^^ - -

EXTENT OF T>€ USEFUL SPECTRUM - 188 ffC

ERROR FUNCTION

DIU » e.488E -2

/«X = a.2SlE -1 rm —e.236E -i

/CAN « 8.681E -€

SIOiftB 8.487E -2

_jl

jJmL 1 iJl L_ 1 . ,1 .II..I II .k. 1 ik u iJ.I

'PI f ^[yr --.ijiqipjf ^ii|' . !|flr '7'

) I

il... .r li yiiiiiiiLUJiikii4iMii4,y

1 1 I^

1

1 1

EXTENT OF THE USEFUL SPECTRUM - 188 »HZ

Figure 4.19 SNR-x 45 dB, SNR-y 46 dB. Deconvolution with moderate regularization,

n fi = 100 MHz, A = 5. |D(n)ldB, d(k), and e(k). Note filtering of

D(n) ,„ is shown in the upper right hand portion

dB

H(n) to obtain D(n) , of the upper figure, vertical scale - 10 dB/div. commencing at 0.0 dB,

horizontal scale, 25 MHz/Div. , commencing at 0.0 MHz.

57
INPUT l*WEFOW

Diu > e.?eeE i

rvix - e.63aE 2

mn »-e.798E 1

HORIZQKTPL SCALE

Diu = 2ee.ee ns

SIGNM. TO NOISE

46.3 06

MfWNITlDE OF T« INPUT WAUEFORM TOURIER TIMNSFORM

10 DB/UIU

HORIZONTAL SCALE

oiu » 25.ee rtc

tHftrvy *l!^ J*wir

^Hiuiii.. 11 ii 1 .1 1 1 II III 1. n

1 111 lUUl llaKiMJBSU

'!.' [7iLii,|J

1 HTf I

1 1 1 1

Figure 4.20 The modified input waveform, x(k) , and the magnitude of its DFT,

I^('^)Ld> after the addition of 35 dB of noise to yield a SNR-x,

dB

40 dB.

58
OUTPUT MAUEFORM

1

i\

V

\Hi^ - -MM j

Diu « e.sesE

rvnn « a.45eE i

niN >-0.sBeE e

HORIZONTN. SOPLE

oiu * 2ee.ee ns

SIOWL TO NOISE

46.2 06

mCNITUDE OF 1>E OUTPUT MAUEFORM FOURIER TIWNSFORM

10 oe/oiu

HORIZONTAL SCALE

DIU « 25.ee MHZ V\b»^J.£}n'^umpiwiirrpi

1 yui 1. 1 Lii . J i:.l'. JV il n

r

1 r

("'

1'

1 1 1

Figure 4.21 The modified output waveform, y(k), and the magnitude of its DFT,

|Y(n) I , after the addition of 6 dB of noise to yield a SNR-y of dB

40 dB.

59
REftL PART OF TVC DECONUOHJTIOH RESULT

DIU » 9.25eE -2

fW< » e.225E -I

MIN «-0.2S0E -2

HORIZONTAL SCALE

OIU > 288.81 HS

FULL SPECTRUM

ERROR FUCTICN

DIU = 0.431E -4

rWX = 8.391E -3

niN =-e.481E -4

MEAN = e.544€ -6

sicn^ e.433E -4

i.

\w l.k \1 • ' 1 1 iL. \^ VifmwiL

'"!f\\ mM

y1 1

iJAiliL JliiUi.b.

pf llT

=t

.v^.. 1

FIXJ. SPECTRUM

Figure 4.22 SNR-x 40 dB, SNR-y 40 dB. Deconvolution without regularization.

|H(n)|^g, h(k), and e(k). H(n) is in the upper right hand

dB

portion of the upper figure, vertical scale -10 dB/Div. conmencinj

at 0.0 dB, horizontal scale 25 MHz/Div. commencing at 0.0 MHz.

60
REAL PART OF T>€ DECONUXUTION RESULT

DIU » e.250E -2

MAX « 0.23aE -1

«IN »-e. 120E -2

HORIZOKTM. SCALE

DIU - 286.61 HS

I

f)^6

\w \I. n

nik 1

1 '1 \"S V\^

%#g nr^ .1,^^ M^g^, k^j>-

EXTENT OF T>€ USEFUL SPECTRUM

ERROR FUMCTIW

74 MHZ

DIU = «.3ai5E -1

MAX « e.issE e

niH =-«.137E 6

nBW « 8.561E -€

SlOVt- 8.341E -1 ^Jk I. .Jul t.Jll .1, n 'ii' T'l'

I 'p'T'iii'fV'r

• MkM\\AM.i\,di^i\iA\bi

"r'l n\f ^ n n| q- f w

EXTENT OF THE USEFUL SPECTRUM - 74 mZ

Figure 4.23 SNR-x 40 dB, SNR-y 40 dB. Deconvolution with strong regularization,

n fi = 74 MHz, A = 6. |D(n)L^, d(k), and e(k). Note filtering

of H(n) to obtain D(n). |D(n)1^3 is shown in the upper right hand

portion of the upper figure, vertical scale -10 dB/Div. coiranencing

at 0.0 dB, horizontal scale 25 MHz/Div. commencing at 0.0 MHz.

61
a^(k)

e(k) pV-AAA—

^^^ 50 fi

a2(k)

J^ Ground-Plane

50 fi u s(k) -• x(k)

Aa^(k)

e(k)g)-^WU50 n -A33 (k)

^

ZI

Ground-Plane

50 '^ s(k) y(k)

Figure 4.24 The time domain antenna range insertion measurement method.

The reference input waveform x(k) is produced by the system

comprised of the 50 ohm source e(k) , the source antenna a^ (k)

,

the receiving antenna a_(k), and the waveform measurement

system s(k). The measurement output voltage y(k) is produced

when a (k) replaces a„(k).

62
JKPvrr fc^OEFORM

DIU » 8-kft«t 3

MAX * e.22SE 4

MiN =-e.?rwF 3

HDRI2DKTAL SCALE

oiu » 2e.ee ks

SIGHAL TO NOISE 1

74.2 oe \;^'

rtAGNITUDE OF 1>€ IhFUT WACEFORM FOURIER TRANSFORM

10 oe^iu

^

HORIZONTAL SCALE

Diu = 258.ee mz \

\

\r\^

1 \ 1 \/V \

(**V|,

H/" ^ i* y

w

(

« V

1 r

Figure 4.25 The input waveform x(k) and the magnitude of its DFT, |x(n)

63
OUTPUT WACtrORM

DIU » 8.2WF 3

MAX =9.16«F 4

MIN =-«.4<WF 3

HORIZOKTAL SCALE

DIU » ^.ee KS

SIO*iL TO NOISE /

67.4 oe r

tlAGNITUOe OF T>€ OUTPUT WAUEFORM FOURIER TRANSFORM

10 DB/OIU

r^

HORIZONTAL SCALE V

DIU » ZM.96 mZ \

\ 1 S -^

u-Ij u. H^WW-w€%tt- ~i ti^^^%

Jl

Figure 4.26 The output waveform y(k) and the magnitude of its DFT, |Y(n)

64
rWO^ITUOE OF THE TRANSFER FUNCTION

HORIZONTAL SCALE

OIU = 258.88 WZ

J i_

1 Ai b ^ 1 t Pt^i^RA r-^^^^^tM^n

' 1

-

Figure 4.27 The magnitude |H(n) I = |Y(n) I t |x(n)

65
RE^L PART OF T»C OECOHUOLUTIOH RESU.T

Diu - e.zeeE -i

nAx e.ieee e

MiN «-e.4eeE -i

HORIZOKTAL SCALE DIU = 28.88 NS

r^ ^^ —

HI W^

EXTEKT OF 1>E USEFUL SPECTRUH 588 MHZ

(a) Filter Parameters: n Q = 500 MHz, A = 20 |D(n) | is shown in upper right hand corner with scales as in Figure 4.27. Error Parameters: |e(k)jmax = 46,

= 43, [|e(k)

'max y(k) ] X 100 = [46/1600]x 100 = 2.9% max

RE^L PART OF THE DECONyOUJTION RESULT

DIU = e.20eE -1

nAX = 8.160E 8

MIN »-e.4e0e -i

HORIZOKTflL SCALE DIU = 28.88 NS

]r -—

^

EXTEhfT CF THE USEFUL SPECTWJH - 308 ffC

(b) Filter Parameters: n ^ = 300 MHz, A = 100, |D(n) | is shown in upper right hand corner with scales as in Figure 4.27. Error parameters: le(k)| = 80, ' 'max 0=7, [|e(k) y(k) ] X 100 = [80/1600]x 100 = 5%.

' 'max

Figure 4.2f

66
[l-R(n)] 4 U-R^)]

1.0

[1-R2(n)]

|y(o)I

u n72

|E(n)| 4

|y(o)

Y(n)|[l-R, (n)]

4221

|Y(n)|[l-R„(n)]

ii^

n, N/2

Figure 4.29 Filtering increases the absolute error |E(n)

67
f(k) DFT[f(k)] = F(n) IDFT[F(n)]

f (k), F (n)

\K,N-l

f (k))

iN\/ . N-l 2N-1 k

G (n)|

k N-l k n-l n

(k)

k

1

,—

N-l

t a^

2k ,2N-1 k

Figure 4-30 Moving from left to right, successive application of the DFT and

IDFT yield periodic functions for |F(n)| and |f(k)|. A cycle of

the resultant g^ (k) may not be the intended f, (k)

b b

68
f(k)

N = 16

. OQ OO OCtO^OOOoOOO

<»-

i M ' I I M I I

< H « t

' ''' ' ' ''' ' '*' ^ 'i

(a)

k N-1

f (k-N)

M = 32

i4»»»»»»4 4 ^^^^f I I M I I I I I I I M <

I

(b)

M-1

f^Ck) = f(k) f (k-N)

M = 32

-H- I I I

I I

I M I I I 1

1 I

' '

' ' ' '

' '» ^»^ » 4f

(c)

M-1

Figure 4-31. (a) The original sequence, f(k), (dark points) is N units in

duration and is constant for k > k . f (k) is extended

by N units (light points).

(b) The extended f (k) is delayed by N units to yield f(k-N).

f(k) and f(k-N) both consist of M = 2N points.

(c) f(k-N) is subtracted from f(k) to yield f2(k).

69
F(n)

|f(o)

N = 16

o o

O O

O 6 H—I—I—I—(—i—4

%—I

\—I—I—I—»- -^ n MfiK- N-l = 15

F^Cm)

2 F(o) --0 6

M = 32

<^<4i4 i 4t4 ' ^i » '4 '^ j ° »'4 ' »'»' (^ i4' »|

-»|i I*- M

2 M^16 M-1 = 31

Figure 4-32. The magnitude of the spectra for F(n) and F (m) The discrete

frequency spacings for F(n) and F (m) are fi and fi/2 due to the

time windows NT and 2NT, respectively. M=2N, thus the total

number of non-zero values in F (m) is the same as in F(n).

70
x(k)

/-

N-1

h(k>

^

y(k)

N-1

. rx^ y 2N-1

h(k)

NI

2N-1

a=l N

4 2

I

\ 2N-]-1

a=l

h(k)

N 2N-1

2

Figure 4-33 Deconvolution using converted step-like waveforms yields _ a pair of impulse responses.

71
KfiL HKT OF 1>€ Dffn»ini.UTION ICSULT

DIU « 9.2SeC -2

MAX e.225E -1

WIN —e.zseE -2

HORIZQKTM. SCALE

DIU > 206.91 HS

1 ..- V

\—

\^^

sH| mk

k

1

\

-I

m

v\ ^-^_^

qsjyw, FUNCTION

DIV > d.4e6E -2

MAX e.l73E -1

riN —e.233E -1

AEAN e.ie4E -5

siQ¥»= e.34ee -z

1 .1 • 1 ill 1

1 1 . 1, . 1 . 1. il

I la^^rri^iil; i-iD^j^ivi.:! hii l:V-Vfl-.rilir:-'-/-^ ; J^ "in fll"! 'fl' 'If

^ n 1 ' rivB

1 T 1

1

1

REQULPRIZATION COEFFICIENT lee.eee

Figure 4.34 SNR-x 75 dB, SNR-y 46 dB. Deconvolution with moderate regularization,

Y 100. |D^(n)|, d^Ck), and e(k). Note filtering of H (n) to obtain

D^(n). |D^(n)|dB is shown in the upper right hand portion of the upper

figure, vertical -10 dB/Div. commencing at 0.0 dB, horizontal 25 MHz/Div.

commencing at 0.0 MHz.

72
VEflL PART OF T>C OECCNUOLUTIGN RESULT

OIU eZSfiE -2

fIftX - e.223E -1

I

\

niN »-e.258E -2 \,> s

H,

ju»a4 .1

HDRIZOKTM. SCALE

1 flL1 OIU - 2ee.ei ns \

-HTpy

v\

^•^^

ERR:5R FUS^TTW

l-TU « A g^pgy -9

(*>* » d.2J^ -1

1 .h 1 ,1. /ijxi »-e.280E -1

fCAH « a.iesE -5

SIQ1A- e.4fl6E -2

PR'H»I[PIHIilWl'^''fl.ir'TlJ

1

:

'—[ —TT

1 • r

1 1 r

REGULAKIZATION CQEFFICIEKT 9Be.eee

Figure 4.35 SNR-x 75 dB, SNR-y 46 dB. Deconvolution with strong regularizatlon,

Y

= 500. [d (n)l, d (k) , e(k). Note filtering of H^(n) to obtain

|d (n)|. |d (n)|dB is shown in the upper right hand portion of the

upper figure, vertical -10 dB/Div. conmencing at 0.0 dB, horizontal

25 MHz/Div. commencing at 0.0 MHz.

73
RE)^ PART GF T>C OECONUOLUTION RESULT

Diy » e.258e -2

MAX « e.225E -1

MiN —e.zseE -2

HORIZOKTM. SCOLE

Diu > zee.ei ns

Diy « e.l4€E -1

rAX » e.733E -1

riN »-e.727E -1

tAH - e.lB9E -5

SI01ft= e.649E -2

^\N \^

V^

V\

\ _%Ll

V\

H •

TM

W)1^

IfM

j4uJ

«^

REOJU^IZftTlON COEfFIClEHT 9999.980

Figure 4.36. SNR-x 75 dB, SNR-y 46 dB. Deconvolutlon with very strong regularization,

Y= 10,000. |D^(n)|, d^Ck), e(k). Note filtering of H^(n) to obtain D (n

|D^(n)|dB is shown in the upper right hand portion of the upper figure,

vertical, -10 dB/div. commencing at 0.0 dB, horizontal, 25 Mllz/div.

commencing at 0.0 MHz.

74
x(k)

>

•Ho

>

6

/ "

1 /

/

/

, / /

50 ps/Div.

y(k)

>

•HQ

?

OCO

^ ^^ -

/

/

J /

/f y\

50 ps/Div.

x(k)

>

•HQ

>

eo

i

£ Lw-.ik«

y(k)

>o

o

en

100 ps/Div. 100 ps/Div.

lX(n)

>

•Ho

U

Y(n)

pqT3 : :::]

- -44-

1 i • i-l-J 4 I* 1 nil 1

5 GHz/Div.

5 GHz/Div.

Figure 4.37 The step-like insertion waveforms x(k) and y(k), and their duration

limited counterparts x(k) and y(k), respectively. |x^(n)| and

|Y (n) I are the magnitudes of the DFT's of x(k) and y(k),

respectively.

75
Y =

Y

= 10

Y

= 10

Y = 10

h

Ml

1

l\

1

' • 1 _L - 1 i il ll il 1

-4k

-1--V4-4^

1B---H- ji ..J.

it

.i, 1 :::::::::::

^^ti.... . ..I...I..

Figure 4.38 The effect of filtering, |d (n) | = |y (n)/X (n)| R(n) . Y = no filtering, Y

= 10 very strong filtering. Vertical scale 10 dB/Div.

horizontal scale, 5 GHz/Div. The regularization filter, R(n) is

real. When Y = 0, |d (n)| = |h (n)|.

76
H(n)

5 GHz/Div.

D^(n)

> _ . ..--,._

•HP

pq 1

1 ti ---[J-- . o llllilllllllllllllllllllllllllll iiiaiiiiiiiiiiiiiiiiiiiiiiiii

5 GHz/Div.

d^Ck)

>Q

>O

E

11 3 2

11 i \lu

[ ^^2^^) r Iw ^ t- tX\ ^ tt

i 4

VD L o J \ ^

z ^:>,^__n

J

100 ps/Div. 50 ps/Div.

e(k)

>

•HQ

> i n^Sjl

t P1 t If'"* H

. _ , ,__

50 ps/Div.

Figure 4.39 The waveforms and spectral magnitudes for very strong

regular izat ion, y = 10 .

77
5. TIME DOMAIN DECONVOLUTION

5.1 Introduction

Time domain deconvolution methods operate directly on the time domain data rather than on transform

domain data. There are certain advantages in time domain operations. Because the physical manifestation

of the signal (pulse, transient, etc.) occurs with the evolution of time, the time domain is the natural

or physical domain of the signal. Direct operations on a physical signal involve operations on real

quantities. The causal and stable nature of a signal is readily recognized in the time domain and can be directly preserved under approximations to, and operations on the signal. Also, convolution and

deconvolution processes implemented in the time domain only employ data corresponding to times such that

k ;£ k for the computation of y(k. ) and h(k^), respectively, in

y(k) = x(k) * h(k) (5.1-1)

h(k) = y(k)(l/*) x(k) (5.1-2)

The presence of noise will ultimately yield undesirable results (probably unstable) in a time domain

deconvolution calculation; this fact is illustrated later in the discussion. However, filtering or smoothing of the time domain data by direct operations on the data can produce stable deconvolution

results, i.e., estimates of h(k). The apparent promise of determining h(k) from the time domain

quantities x(k) and y(k) by direct operations in the time domain is a sufficient stimulus to try and possibly not fail.

The time domain methods presented here are not as easy to apply in comparison to the frequency

domain methods presented in Chapter 4, nor are the results as smooth. In the authors' view, the time

domain methods are still in the experimental stages of development while the frequency domain methods

can be used with a significantly greater degree of confidence on the part of the user. However, the promise of the time domain method remains to spur further development studies.

In the remaining sections of this chapter the following topics will be discussed:

1. Classical time domain deconvolution

2. Iterative time domain deconvolution

3. Regularization or filtering in the time domain

4. Experimental results of classical deconvolution with regularization applied to the

RG 58 C/U transmission line problem of Chapter 4.

5.2 Discrete Classical Time Domain Deconvolution

5.2.1 Convolution

The discrete classical time domain deconvolution method is an algebraic process which is the inverse

of the convolution algebraic process [26]. For the moment consider the multiplication of the two polynomials f . (x) and f.(x).

f,(u) = a + a,u + a„u^ + ••• + a^, ^u^"-*- (5.2.1-1)

1 o 1 2 N-1

f„(u) = b + b,u + b^u^ + ••• + b„ .u^"-*" (5.2.1-2)

I O 1 / N-1
The expression for the coefficients of the product

f3(x) = f^(x) f^Cx) (5.2.1-3)

is identical to the convolution of the N term sequences

h(k) = a^, a^, a^, ••• , a^_^ (5.2.1-4)

x(k) = b^, b^, b^, ••• , h^_^ (5.2.1-5)

and is given by

f y(k) = h(k) * x(k) = Coefficients of the product, f (u)f„(u)

k-1

^ i?0

""^^^"^ (5.2.1-6)

In the notation used in the earlier chapters, the convolution result for y(k) is given by

k-1

y(k) = ^ h(i)x(k-i) (5.2.1-7)

i=0

Because convolution corresponds to multiplication of two polynomials, deconvolution corresponds to division of two polynomials. These views are clearly seen in the Z-transform representation which is discussed in section 5.2.2. For the causal signals h(k) and x(k), both are equal to zero for k < 0. Expansion of (5.2-7) for causal signals leads to a result in which the k = term is zero; hence the resultant sequence for y(k) commences with k = 1, figure 5-la. Accordingly, deleting the initial zero values h(0) and x(0),

(5.2-7) may be written as

k

y(k) = X^hCi) x(k+l-i) (5.2.1-8)

i=l

which yields the set of equations for the y sequence values.

yd) = h(l) x(l) (5.2.1-9)

y(2) = h(l) x(l) + h(2) x(l) (5.2.1-10)

y(3) = h(l) x(3) + h(2) x(2) + h(3) x(l) (5.2.1-11)

and so on. Notice that y(k ) is computed from earlier values of h(k) and x(k)

i.e., k £ k .

5.2.2 Z-transform Representation

Some of the properties of convolution and deconvolution are more easily developed in the Z-transform

domain [29]. Before proceeding with the discussion, some of the properties of the Z-transform will be briefly described. Given a pair of consistent (exact) discrete waveforms x(k) and y(k), each one of them

can be transformed to the complex z domain by the Z-transf ormation

79
F(z) = Ef(k) z " (5.2.2-1) k=0

that is, f(k) is transformed to F(z); hence

CO X(z) = 2] x(k) z"^ (5.2.2-2)

k=0

OO Y(z) = XI y(k) z'^ (5.2.2-3)

k=0

Furthermore, the discrete convolution

y(k) = x(k) * h(k)

transforms to Y(z) = X(z) H(z) (5.2.2-4)

The ratio of the two polynomials for Y(z) and X(z) yields H(z)

,

Z y(k) z ^

„, , k=0

H(z) = (5.2.2-5)

E x(k) z~^

k=0

which is another polynomial in z,

H(z) = ^ h(k) z ^ (5.2.2-6)

k=0

Consequently, as it was pointed out in Section 5.2.1 the coefficients of the convolution and deconvolution sequences are generated according to the rule for polynomial products and quotients,

respectively.

5.2.3 Deconvolution

A set of equations for generating the h sequence can be obtained from eqs. (5.2.1-9) through

(5.2.1-11) (and so on) by first rearranging them to

h(l) x(l) = yd)

h(2) x(l) = y(2) - h(l) x(l)

h(3) x(l) = y(3) - h(l) x(3) - h(2) x(2) and so on. In summation form,

k-1

h(k) x(l) = y(k) - 22 h^i) ^^^ + 1-i^ (5.2.3-1)

i=l

Dividing x(k) and y(k) by x(l) gives

x(k) = jg^ (5.2.3-2)

80
y(k) = ^

(5.2.3-3)

where

x(l) i (5.2.3-4)

x(l) = 1 (5.2.3-5)

Consequently, h(k) can be written as

_ k-1

h(k) = y(k) - V h(i) x(k+l-i) (5.2.3-6)

1=1

Briefly, (5.2.3-6) generates the following sequence for h(k);

h(l) = 7(1) (5.2.3-7)

h(2) = 7(2) - h(l) x(2) (5.2.3-8)

h(3) = 7(3) - h(l) x(3) - h(2) x(2) (5.2.3-9)

h(4) = 7(^) - h(l) x(4) - h(2) x(3) - h(3) x(2) (5.2.3-10)

h(5) = 7(5) - h(l) x(5) - h(2) x(4) - h(3) x(3) - h(4) x(2)

(5.2.3-11)

Inspection of the values of h(k) on a line to line basis shows that the earlier values of the impulse response sequence, e.g., h(l), h(2), ••', h(k -1), are used to calculate the specific value h(k ). Because the h-values can be expressed in terms of the input and output sequence values, x and y, the h values are actually being computed from the x and y values. Note that h(l) can be eliminated from the expression for h(2) by using the first relation for h(l) and h(2).

h(2) = 7(2) - [7(1)] x(2)

Similarly, for h(3) there results,

h(3) = 7(3) - 7(1) x(3) - [7(2) - 7(1) x(2)] x(2)

continuing the process obtains the partial set of equations for h(k)

,

h(l) = 7(1) (5.2.3-12)

h(2) = [7(2) - 7(1) x(2)] (5.2.3-13)

h(3) = [7(3) -7(1) x(3)] - [7(2) -7(1) x(2)] x(12) (5.2.3-14)

h(4) = [7(4) - 7(1) x(4)] - [7(2) - 7(1) x(2)] x(3) - [7(3) - 7(1) 7(3)] x(2) + [y(2) - yd) x(2)] x^(2) (5.2.3-15)

h(5) = [7(5) - 7(1) x(5)] - [7(2) - 7(1) x(2)] x(4)

- [7(3) - 7(1) x(3)[ x(3) + [7(2) - 7(1) x(2)] x(2) x(3)

- [7(4) - 7(1) x(4)] x(2) + [7(2) - 7(1) x(2)] x(3) x(2)

+ [7(3) - 7(1) x(3)] x^(2) - [7(2) - 7(1) 7(2)] 7^(2) (5.2.3-16)

81
h(l) ==

h(2) == a(l) A(2)

h(3) == a(l) A(3) + a(3) A(2)

and so on. Defining the quantities,

A(i) = y(i) - yd) x(i) (5.2.3-17)

aCl) = 1 (5.2.3-18) r-1 _

ct(r) = - 2 "(i) x(r+l-i) (5.2.3-19)

i=l

enable h(k) to be written as

(5.2.3-20)

(5.2.3-21)

(5.2.3-22)

h(4) = a(l) A(4) + a(2) A(3) + a(3) A(2) (5.2.3-23)

h(5) = a(l) A(5) + a(2) A(4) + a(3) A(3) + a(4) A(2) (5.2.3-24)

and so on. Explicitly, in terms of (5.2.3-16) through (5.2.3-19), h(k) is given by

k-1

h(k) = yd) + 2 a(i) A(k+l-i) (5.2.3-25)

i=l

k-1 • _ _ _

= yd) + 22 "*^^^ [y(k+l-i) - yd) x(k+l-i)] (5.2.3-26)

i=l

which is the general result for h(k) in terms of x(k) and y(k). To gain further insight into this result for h(k), (5.2.3-26), consider the signals shown in

figure 5.1. x(k) and h(k) are given signals while y(k) is calculated from x(k) * h(k), figure 5.1a. Subsequently, the deconvolution result, h (k) , is calculated from y (k) (l/*)x(k) . Note that both y(k) and

h , (k) commence at k = 1 while x(k) and h(k) commence at k = 0. From (5.2.3-17) through (5.2.3-19) using

the data of figure 5.1, there results

A(2) = 7(2) (5.2.3-27)

(5.2.3-28)

(5.2.2-39)

(5.2.3-30)

and so on. Next,

(5.2.3-31)

(5.2.3-32)

(5.2.3-33)

-x(3) +xh)] x(2) = 2x(2) - x-^(2) (5.2.3-34)

and so on. Using the above values, calculate the deconvolution of y(k) and x(k) to obtain h.(k). To

d insure clarity, the deconvolution result for h(k) is denoted by hj(k). The terms in h , (k) will be

d d uniquely related to those of the given sequence, h(k). From (5.2.3-20) through (5.2.3-24),

k-1

^d'^^^ " E a(i)A(k+l-i) (5.2.3-35)

i=l

The first four terms are

82

A(3) = y(3)

A(4) = 7(4)

A(5) = y(5)

a(l) = 1 a(2) = x(2)

a(3) = -x(3) + x2(2)

a(4) = x(2) x(3) -
h^(l) ==

h^(2) == h(l)

h^(3) == h(2)

h^(4) == h(3)

hjjd) = (5.2.3-36)

h^(2) = a(l)A(2) = y(2) (5.2.3-37)

h^(3) = a(l)A(3) + cx(2)A(2) =7(3) - x(2) y(2) (5.2.3-38)

h^(4) = a(l)A(4) + a(2)A(3) + a(3)A(2)

= 7(A) - x(2) 7(3) + [-x(3) + x^(2)] 7(2) (5.2.3-39)

Upon substituting the values of y(k) from the convolution (5.2.1-8) into (5.2.3-36) through (5.2.3-39)

and recalling that x(k) = x(k)/x(l) and y(k) = y(k)/x(l), the following values for h (k) are obtained.

(5.2.3-40)

(5.2.3-41)

(5.2.3-42)

(5.2.3-43)

Inspection of these results shows that h (k) is equal to h(k-l), figure 5.1b. Figure 5.2 shows the result of classical deconvolution using two noise free signals, x(k) and y(k).

In this case x(k) and h(k) were specified as

x(k) = e

"^ (l-e~'^ ); k = 0,1, •••,49 (5.2.3-44)

h(k) = e'-^^^d-e'"'^) ; k = 0,1, •••,49 (5.2.3-45)

which, in turn, specified y(k) as

y(k) = x(k) * h(k) k = 0,1, •••,49. (5.2.3-46)

The data for y(k) was calculated using (5.2.1-8). The resultant y(k) was used with x(k), (5.2.3-44),

in (5.2.3-26) to calculate the deconvolution.

h(k) = y(k)(l/*) x(k) k = 0,1, •••,49 (5.2.3-47)

In applying (5.2.3-6) or (5.2.3-26) it should be kept in mind that for k equal to zero, h(k) is always

zero (causal system); consequently, the first term, zero, must be added to the sequence for h(k). In the present example, (5.2.3-6) or (5.2.3-26) yields a sequence of 49 terms, and including zero as the first term in the sequence gives a total of 50 terms in the sequence.

5.2.4 The Effect of Errors in Classical Deconvolution

The effect of errors or noise on the deconvolution process can be seen by considering (5.2.3-6) and

its partial expansion (5.2.3-7) through (5.2.3-11) in the presence of errors. For the moment, recall

the deconvolution example given above in Section 5.2.3. The noise free signals x(k) and h(k) , (5.2.3-44)

and (5.2.3-45), respectively, were used to calculate a noise free y(k), (5.2.3-46); this y(k) was con- sistent with x(k) and h(k) through the convolution equation (5.2.1-8). Now imagine that one value of the y(k) sequence is changed to an incorrect value; let y(l) be altered to y(l),

a.7(1) = 7(1) + e (1) (5.2.4-1)

83
where e (1) is the error added to the correct yalue, y(l). Putting y(l) into (5.2.2-7) through (5.2.2-11)

yields a partial set of the erroneous h(k) s, h(k),

h(l) = h(l) + e (1) (5.2.4-2)

h(2) = h(2) + e (1) [-x(2)] (5.2.4-3)

h(3) = h(3) +e (1) [-x(3) +x^(2)] (5.2.4-4)

h(4) = h(4) + e (1) [x(4) - 2x(2) x(3) - x-^(2)] (5.2.4-5)

h(5) = h(5) +e (1) [x(5) +x^(3) +x^(2) +x^(2)] (5.2.4-6)

and so on. The effect of the error, e (1), is increased with increasing k; the error in h(k) is simply

e^(k) = h(k) - h(k) (5.2.4-7)

h

Here, the set of errors, e (k) , in h(k) is due to the single error in y(k) , e (1). In practice, the data values for x(k) and y (k) will each be corrupted by noise; consequently, the errors in h(k) will be

dependent upon the sets of error for both x(k) and y(k). Consequently, one car be very sensitive to the errors and most probably will be an unstable sequence.

The analysis can be extended to a set of errors in be free. Let the erroneous values of y(k) be given by y(k).

dependent upon the sets of error for both x(k) and y(k). Consequently, one can imagine that h(k) will

;ry sensitive to the errors and most probably will be an unstable sequence.

The analysis can be extended to a set of errors in both y(k) and x(k) . First, let x(k) be error

y(k) = y(k) + ey(k) (5.2.4-8)

Putting y(k) into the general result for h[x(k), y(k)], (5.2.3-26) obtains

h(k) = h(k) + e (1) + y :t(i) [e (k+l-i) - e (1) x(k+l-i), (5.2.4-9)

Vy

*.- V V i^i y y

The error in h(k) is e (k) and is given by

e^(k) = h(k) - h(k) (5.2.4-10)

h k-1 _

e,(l) + "Z, "(i) [e (k+l-i) - e (1) x(k+l-i)] (5.2.4-11)

y i=l

Next, include the errors in x(k) by replacing x(k) and a(k) by

x(k) = x(k) + e (k) (5.2.4-12)

X

%

a(k) = a(k) + e^(k) (5.2.4-13)

The final result for the erroneous h values is given by

_ k-1 _ _ _ k-1

h(k) = yd) + J2 '-*(i) [y(k+l-i) - yd) x(k+l-i)] - 2^ "(i) yd) e^(k+l-i)

i=l i=l

k-1 _ y e (i){y(k+l-i) - yd) [x(k+l-i) + e (k+l-i)]) + (con't, next page)

i= 1

''

84
k-1 _

+ e (k) + X] [a(i) + e (l)]{e (k+l-i)-e (1) [x(k+l-i) + e (k+l-i) ] } (5.2.4-14)

q y

The error in h(k) is then

e^(k) = h(k) - h(k) (5.2.4-15)

h k-1 _ k-1 _ _ _

e,(k) - Y ^(i) y(l) e (k+l-i) + Y e„(i) {y (k+l=i)-y (1) [x(k+l-i)+e (1+1-i) ]

}

y i=l i=l

K-1 _

+ Yl [•(i)+e (i)] {e (k+l-i)-e (1) [x(k+l-i)+e (k+l-i)]} (5.2.4-16)

i=i 1 y y ^

To illustrate the effect of noise on the deconvolution process, noise was added to the noise-free

signals (5.2.3-44) and (5.2.3-46) so that each signal possessed a 40 dB signal to noise ratio (as defined in Section 4.1); the result is h(k) , (5.2.4-14), and is shown in figure 5.3. It is clear that the errors (noise) lead to an unstable result which bears no resemblance to h(k), figure 5.2. Similar

results are also obtained for the iterative deconvolution method which will be discussed in section 5.3

5.2.5 The Critical Dependence of Classical Deconvolution

on the First Point in x(k).

In the practical application of the classical time domain method to the deconvolution problem, the choice of the first term in the input sequence, x(l), is critical. If x(l) is ill-chosen the resultant

deconvolved result for h(k) will not converge. To gain an understanding of the origin of this feature,

consider the Z-transforms of the two known sequences x(m) and y(n) having their durations denoted as

M' and N' , respectively. Thus,

X(z) = x(l) z"-"- + x(2)z~V •••+ x(M')

= z~'''[x(l) z^* + x(2) z^'"-^+ •••+ x(M+l)]; M'-1 = M (5.2.5-1)

Y(z) = yd) z"-^ + y(2) z"^ + ••• + y(N')

= z~^[y(l) z^ + y(2) z^"""" + ••• + y(N+l)] N'-l = N (5.2.5-2)

and

respectively, where

N>M (5.2.5-3)

The Z-transform of the impulse response corresponding to x(k) and y(k) is

nil) = Y(n)/X(m) (5.2.5-4)

When N equals M, H("".) becomes H(n) and is a constant.

H(n) = ^^ = constant =C; N = M (5.2.5.5)

A(,n)

The inverse Z-transf orm of H(n) is then the impulse sequence C, 0, 0, ••• in which the first term is C and is followed by N-zeroes.

By synthetic division H(z) is given by

85
H(z) = -^!^z^-.ry(2)-^^^%^1^

N-M-1 N-M x(l) \ y^^'

5^(2^) J x(l)

z ' L

i^^n^ _ y(i)x(3) 1 _ r . _ y(i)x(2) 1 2L(2i j z^-^-^ .. j

j[y^^^ x(l) J [y(2) ^^^^ J ^(^) j ^(^) + ^ (5.2.5-6)

Using the variables x(k) and y(k) which are the x(k) and y(k) normalized to x(l), H(z) can be written

as

H(z) = 7(1) + [7(2)-7^(l)x(2)]z"^ + {7(3)-7(l)x(3)-[7(2)-7(l)x(2)]x(2)}z"^ + •

(5.2.5-7) Now consider that the point x(l) is corrupted by an error, e , and thus changed to x(l)

x(l) = x(l) + e (5.2.5-8)

while all other values of x(k) and y(k) are correct. The normalized x(m) and y(n) would be given by

-/ ^ 1 x(2) _. x(M')

x(.m) = 1, —7TT7 > » —7T^C~Z x(l)+ e x(l) + e X X

= 1, x(2) + e (2), •••, x(M') + e (M') (5.2.5-9)

-( ^ = y(l) y(2) ... y(N')

^^ ' x(l) + eXX' x(l) + e ' ' x(l) + eX

= yd) + e (1), y(2) + e (2), •••, y(N') + e (N') (5.2.5-10)

respectively. Then the only correct term in x(m) would be x(l) which equals x(l) = 1. The y(n)

sequence would be scaled to the erroneous value x(l) + e . Consequently, the ratio for H(z) would also

be erroneous,

J^(z) = li^ = 7(1) + [7(2) - 7(1) x(2)] z " +

{ y(3) - yd) x(3)

X(z)

- [7(2) - yd) x(2)] x(2)} z +••• (5.2.5-11)

The inverse Z-transform of H(z) gives the erroneous sequence h(k), (5.2.4-14), and an error in h(k) of

, (5.2.4-16).

Consequently, it has been shown that the error terms, i.e.,

e^(k) : e^(2), e^(3) , •••

-a^^)

= ^a(2),

e^(3), •••

e (k) : e (1), e (2), ••• y y y

86

e^(k), (5.2.4-16)

h
e^(k) : e^(2), e^(3), •.•

have been generated by the single error in x(l) through the normalization processes (5.2.5-9) and (5.2.5-10). Thus the selection of x(l) is critical to the deconvolution process.

5.3 Discrete Iterative Time Domain Deconvolution

5.3.1 Iterative Method

By using an iteration process it is also possible to solve the deconvolution equation by an iterative convolution process [4,5,27,28] proposed by P. H. van Cittert. [28]. The principle of this

method is to form successive approximations, »h.(k), to the unknown system impulse response, h(k) , using

the convolution equation,

y(k) = h(k) * x(k) (5.3.1-1)

For physically realizable signals the duration of x(k) is less than y(k). As a first approximation

to h(k), h^(k), use y(k),

h^(k) = y(k) (5.3.1-2)

Using h^ (k) in the convolution equation (5.3.1-1) gives the result y (k)

y^(k) = h^(k) * x(k) (5.3.1-3)

y. (k) differs from the true result y(k) by the error, [Ay(k)]

[Ay(k)]^ = y(k) - y^(k) (5.3.1-4)

which corresponds to an error in h(k) of

[Ah(k)]^ = [Ay(k)]^ (1/*) x(k) (5.3.1-5)

correcting h (k) by adding [Ah(k)] yields a second approximation to h(k), h„(k).

h2(k) = h^(k) + [Ah(k)]^ (5.3.1-6)

= h^(k) + [y(k) - h^(k) * x(k)] (5.3.1-7)

Consequently, the i-th approximation is related to (i-l)th approximation by

h.(k) = h^_^(k) + [y(k) - h^_^(k) * x(k)] (5.3.1-8)

When ^

h.(k) * x(k) = y(k), (5.3.1-9)

1 then

h^(k) = h^^^(k) = h(k) (5.3.1-10)
Here, the y(k) and x(k) are considered to be noise free. Of course, they are not in practice, and the noise does limit the direct application of the iterative method.

Continuing the discussion on the iterative method, apply the Z-transformation to (5.3.1-8).

The result is

Y(z) - H^(z) X(3) + H^(z) = H._|_^(z) (5.3.1-10)

The form of (5.3.2-7) parallels that of (5.3.1-8) with the convolution operation replaced by a product.

The iteration process starts with the first estimate for H(z) being chosen as Y(z),

H^(z) = Y(z) (5.3.1-11)

Putting H (z) into (5.3.2-4) yields for U^(z)

n^(z) = Y(z)[2-X(z)] (5.3.1-12)

Repeating the process for the next values yields

n^iz) = Y(z) + H2(z)[l - X(z)] = Y(z) {1 + [2 - X(z)][l - x(z)]}

(5.3.1-13)

H,(z) = Y(z) {1 + [l-X(z)] + [2-X(z)][l-Y(z)]^} (5.3.1-14)

4

H^(z) = Y(z) {1 + [l-X(z)] + [l-X(z)]^ + [2-X(z)][l-X(z)]^} (5.3.1-15)

H,(z) = Y(z) {1 + [l-X(z)] + [l-X(z)]^ + [l-X(z)]^

6 + [2-X(z)][l-X(z)]^} (5.3.1-16)

and so on. The result for H.(z) is a power series of the form

[± + B^ X(z) + B2 X(z)^ + ••• + B._^(z)'-

-^J H^(z) = Y(z) 1 i + B, X(z) + B^ X(z) + ••• + B._, (z) \ (5.3.1-17)

whose inverse Z-transform is

h.(k) = iy(k) + B^y(k) * X(k) + B2y(k) * X(k) * X(k) + •••

+ B. , y(k) '•= x(k) * x(k) " •••x(k) (5.3.1-18)

1-1 "

Hence, it is seen that the direct discrete iterative method requires a large number of convolution

operations which in turn will require a large amount of computation time.

Using the x(k) and y(k) of (5.2.3-44, 5.2.3-46), the direct discrete iterative deconvolution method

was employed to calculate h(k); the results for 500 and 1000 iterations are shown in figures 5.4 and 5.5 where it is seen that the latter result is identical to the result obtained by the direct discrete

classical deconvolution method. However, the discrete iterative method required 45 minutes of computation time as compared to 10 seconds for the classical method.
5.3.2 The Effect of Errors in Iterative Deconvolution

If y(k) contains an error, say for k = 0, then the erroneous output sequence y(k) is

y(k) = 6^(0) +y(0), yd), y(2), ••• (5.3.2-1)

where e (o) is the error in y(0). The Z-transform of y(k) is

Y(z) = e (0) + y(0) + yCDz"""" + y{2) z"^+ ••• (5.3.2-2)

= e (0) + J2 y(k) z"^ (5.3.2-3)

^ k=0

= e (0) + Y(z) (5.3.2-4)

Using the erroneous output Y(z) in (5.3.1-17) yields the erroneous i-th approximation to the deconvolution

result

H^(z) = [e (0) + Y(z)] [i + B^X(z) + B^X^ {z) + ••• + B ._^x'^"-'-(z) ] (5.3.2-5)

= + H.(z) +e (0) F.[X(z)] (5.3.2-6)

where the function of X(z) is defined by

2 ^-1 F^[X(z)] = [i + B^X(z) + B2X (z) + ••• + B^_^X(z)] (5.3.2-7)

Consequently, the effect of the error e (0) increases with increasing k. This result is similar to the error result for the classical deconvolution technique with a single error (5.2.4-6).

For a set of errors, e (k) , in y(k), i.e., for an erroneous output sequence of the form

y(k) = e (k) + y(k) (5.3.2-8)

the Z-transform of the sequence is

Y(z) = E (z) + Y(z) (5.3.2-9)

Using (5.3.2-9) in (5.3.1-17) yields the erroneous i-th approximation to the deconvolution result

H^(z) = [E (z) + Y(z)] F^[X(z)] (5.3.2-10)

= H.(z) + E (z) F.[X(z)] (5.3.2-11)

whose inverse Z-transform yields the erroneous time domain sequence, h.(k),

h.(k) = h.(k) + e (k) * f.(k) (5.3.2-12)

where f . (k) is the inverse Z-transform of F.[X(z)]. h.(k) is the y-noise or error corrupted result for the noise free i-th iteration approximation, h.(k), to the system impulse response, h(k).

If x-noise or errors, e (k) , are also present they transform to E (z) ; and then X(z) becomes X(z)

,

89
X(z) = X(z) + E^(z) (5.3.2-13)

and the corrupted i-th iteration, is given by

H.(z) = H.(z) + Ey(z) G.[X(z) + E^(z)] (5.3.2-14)

G.[X(z) + E (z)] is shown in expanded form in (5.3.2-16). The error in H.(z) is then

x^ X X

Ej^^(z) = H^(z) - H.(z)

= E (z)G [X(z)+E (z)] (5.3.2-15)

y 1 X

= E (z) ] i + B [X(z)+E (z)] + B„[X(z)+E^(z)]V •••+B. , [X(z)+E (z)]^"^}

y ' -L X z X i—x X (5.3.2-16)

In the time domain, the i-th corrupted iteration, h.(k), for h(k) is given by the inverse Z-transform

of (5.3.2-14).

h^(k) = h^(k) + Ey(k) * g^[x(k), e^(k)] (5.3.2-17)

where g.(k) is the inverse Z-transform of G.[X(z) + E (z)].

Figure 5.5 shows the result for h.(k) with i=750 (iterations) using noise corrupted values of

•^ """% (5.2.3-44) and (5.2.3-46) for x(k) and y(k), respectively,

x(k) = x(k) + e (k) (5.3.2-18)

y(k) = y(k) + e (k) (5.3.2-19)

where e (k) and e (k) were such that the SNR of x(k) and SNR of y(k) were both 40 dB. Note that the

x y result is unstable and bears no resemblance to the system response h(k), (5.2.3-45). Figure 5.5 is very similar to figure 5.2, the result for noise corrupted classical deconvolution under the same x and y signal to noise ratios (40 dB)

5.4 Time Domain Regularization

5.4.1 Averaging Filter

In sections 5.2 and 5.3 two direct methods of deconvolution were presented. By direct it is meant

that the mathematical method is directly applied to the data x(k) and y(k). As we demonstrated in sections 5.2 and 5.3 the presence of noise or errors in the data destroys the consistency between x(k)

and y(k); that is to say, the computed data x(k) and y(k) are no longer mathematically related to each

other through the impulse response of the system, h(k) . Thus, the direct application of the mathematical

methods to the computed data x(k) and y(k) leads to an erroneous and divergent or unstable result, h(k)

90
Both the classical and iterative methods applied to noisy signals lead to the same type of Incorrect

and divergent results.

By filtering the noisy waveforms x(k) and y(k) both in the same prescribed way, the functions x(k)

and y(k) can be smoothed or regularized so that when the mathematical deconvolution method is applied to them, the process yields a stable convergent result, d(k), which approximates h(k). The smoothing is accomplished by a regularization filter. Figures 5.7 and 5.8 show the resultant estimates for h(k)

obtained by first filtering the corrupted signals x(k) and y(k) which were used in earlier examples.

Figure 5.7 was obtained by using the filtered signals in the discrete classical deconvolution method,

(5.2.3-6) or (5.2.3-26). Simarily, figure 5.8 was produced by the regularized discrete iterative deconvolution method, i.e., filtered x(k) and y(k) in equation (5.3.1-8). Comparison of figures 5.7 and 5.8 shows that the two methods yield very similar results, but again, the computation time for the iterative

method was 45 minutes compared to 10 seconds , for the classical deconvolution method. The notation d(k)

denotes an estimate of the true impulse response h(k) which was estimated from the filtered or regularized

data x(k) and y(k).

^\, 'X, The regularizzation filter applied to x(k) and y(k) was a simple one, a four-point averaging filter: the value of, say x(k) , at the point k was set equal to the average or mean value of the values of x(k) , x(k+l), x(k+2), and x(k+3); i.e., x(k) was set equal to x(k)

.

N ^(k) = ^ E ^(J) (5.4.1-1)

j=k

where N = k + 3. Without the filter (5.4.1-1), the S/N ratios for x and y of 40 dB led to failure of the direct deconvolution methods. With the filter (5.4.1-1), the S/N ratios of 20 dB led to failure;

consequently, it could be said that there was a 20 dB enhancement of the deconvolution process. A major task in time domain deconvolution is the synthesis of a suitable regularization filter.

The example described above employed the simple averaging filter (5.4.1-1); however, to achieve improved

estimates of the impulse response, d(k), more sophisticated regularization filters need to be synthesized.

In the next section a more sophisticated and flexible filter is described.

5.4.2 Composite Filtering

5.4.2.1 Introduction

In the last section an averaging filter was briefly discussed and applied to two deconvolution

problems. The averaging filter is equivalent to a system having a rectangular impulse response, h (k)

.

h (k) 1' 1k <k^ (5.4.2.1-1)

0, k < k < N-1

In other words, the impulse response h (k) consists of N discrete values, k^ + 1 of which are equal to

a 1 unity, while the remaining values are zero.

If a given waveform sequence x(k) is convolved with h (k) using (5.2.1-7) or (5.2.1-8), the result

is the same as that produced by the averaging operation, (5.4.1-1).

Presented below is a filter consisting of three operations which are (1) averaging, (2) least-square

first-order polynomial fitting, and (3) Gaussian weighting or smoothing. Each one of these operations is represented by an adjustable parameter. Consequently, the overall composite filter possesses three

91
parameters, L, S, and A corresponding to averaging, fitting, and weighting, respectively.

5.4.2.2 Filter Description

The main characteristic of this filter is that it is possible to adjust its transfer characteristic

over a relatively wide range by adjusting the three parameters L, S, and A. The transfer properties will

be considered in the next section while here the filter structure will be presented. L and S are integers

while A is a positive real number. Let W(k) be the discrete waveform to be filtered.

L is the running averaging parameter and is equal to the time window (or number of points) of the averaging operation. Each value of W(k) , say W(k ), is replaced by the mean value of W(k. ) and the L-1 values of W(k) following W(k ) , i.e.,

New W(k) = L S W(k+i) (5.4.2.2-1)

S is the polynomial fitting parameter and is equal to the time window (or number of points) over which the least square, first order polynomial fitting [30] is performed. The window is shifted along

the waveform in such a way that each point of the waveform W(k) , say W(k. )

, comes S times in the time

window. For each position of the window the S selected points are replaced by a new set of S values

deduced from the fitting. The first time the point W(k ) appears in the window it is transformed to a new value W(k ,1); the second time the value W(k ,2) is obtained and so on until the value W(k ,S) is obtained.

A is the characteristic constant of a Gaussian weighting function which is applied to the S values

of W(k , i)

S New W(k) S W ^(k) = ^ a(i) W(k,i) (5.4.2.2-2)

' i=l

where

i=l

The expanded form of (5.4.2.2-2) is then

_.2, 2 a(i) = a e

" ^ (5.4.2.2-3)

S^ a(i) = 1 (5.4.2.2-4)

-2 -2 2 -2 Wg ^(k) = a I e"'^ W(k,l) + e""*^ W(k, 2) + • • • + e"^ "^ Te ^ W( W(k,S)| (5.4.2.2-5)

= le +e +'"+e

where

-2 ,.-2 o2.-2

(5.4.2.2-6)

Summarizing, the composite filter employs three filtering or smoothing techniques: averaging, least

squares-first order fitting, and Gaussian weighting, with the corresponding adjustable filter parameters

L, S, and A, respectively. The filter is divided into two filters, the L-averaging filter and the S, A- fitting/weighting filter.

92
5. A. 2. 3 Frequency Domain Analysis of the Composite Filter

Frequency domain analysis of the composite time domain filter (or any other time domain filter) is not necessary for the application of time domain filtering. It is presented here to provide the reader

with additional insight to the filtering processes by establishing a connection to the frequency domain

and thus to some conventional frequency domain concepts and terminology.

Upon applying each filter, (5.4.2.2-1) and (5.4.2.2-2), to an impulse (a single data point) the impulse response of each filter can be obtained. The resultant impulse response sequences can then be Fourier transformed to the frequency domain.

The transfer function of the averaging (L-parameter) filter (5.4.2.2-1) is characterized by the cut- off frequency fco whose value decreases as L increases. For example, for N = 1024 points and a total

time window of 200 ns the corresponding values for L and fco are as follows: L= 2, 3, 4, 5, 6, 7, 8, 9 and 10, and fco (MHz) = 500, 340, 255, 205, 170, 145, 130, 115 and 105, respectively.

The transfer function for an L-f liter with L = 3 is shown in figure 5.9A. Also shown is the transfer function for an S, A filter with S = 5, A = 1 and the cascade (logarithmic sum) of the two filters. Note that the L filter cut-off frequency was chosen to occur at the first maximum of the S, A

filter to yield a sharper cut-off characteristic and increased attenuation above fco. The transfer properties of the S, A filter depends upon the two parameters S and A. Figure 5.10

shows how the value of A affects the transfer properties with S being held constant and equal to 5. Similar dependencies on A appear in figures 5.11 and 5.12; however, figure 5.10 provides a systematic

view of the effects of A. In figure 5.10, A is varied over the range of 0.5 through 5.0. For S = 5, the value of A = 2.05

is of special significance. As A is increased from 0.50, the first two minima shift towards each other

and merge together when A equals 2.05. The resultant cut-off characteristic or minimum is very deep,

about 100 dB. As A is increased above 2.05, the cut-off frequency remains the same while the overall

transfer function broadens out and becomes smoother.

In figures 5.11 and 5.12 examples of composite filtering are shown. Figure 5.13 illustrates the

S, A filter cut-off frequency dependence on the values of S and A, and can be used as a design tool for selecting S and A to achieve a given cut-off frequency. The cut-off frequency, fco, can be scaled

with Af; in the present case Af = 50 MHz/div.

In figure 5.13 for a given fco there are multiple choices for the pair of S, A values. For example,

if fco is desired to be 120 MHz, figure 5.13 indicates that two possible choices for S, A are 12, 4.4 and 13, 5.61. Some contours of limited extent are shown for A equal to 1.0 through 6.0 along with

single points for A over the range of 2.51 through 6.49. For 1.0 ^ A 2l 2.5, the cut-off frequencies of the S, A filters in figures 5.9 and 5.10 can be

compared to those in figure 5.13. Also by using a time window of 200 ns (Af = 50 MHz) in figures 5.11 and 5.12, the cut-off frequencies in figures 5.11 and 5.12 can be compared to those in figure 5.13. For example, from figure 5.11 fco = 250 MHz occurs for S, A equal to 6, 2.40. Inspection of figure 5.13 shows that for S = 6, A would be between 2 and 3.

93
5. A. 3 Application of Time Domain Filtering

5.4.3.1 Practical Procedure

Time domain deconvolution is a progressive process (one point after another) . Knowledge or selection

of the first point in the input sequence, x(l), is critical as was shown in section 5.2.5. The selection

o'f x(l) is accomplished by trial and error using the deconvolution algorithm; however, this is not as random an exercise as it may appear to be. The search for x(l) is performed using only a small part of the available data points, say 50 points out of 1024. Typically, a choice is made for x(l) and the deconvolution algorithm is run for 50 points or so. If the result is stable, then the algorithm is applied to increasingly larger numbers of data points until the entire data span is Included, with x(l)

being adjusted until a value is found which provides an overall stable deconvolution result. There are two factors that strongly interact with the selection of x(l) : (1) data offset error and

(2) noise in the x and y data. The offset error correction can be applied before commencing the deconvolution operation (refer to section 4.6.3). Noise in the data can only be eliminated by filtering;

consequently, some combination of offset correction, x(l) selection, and filtering is necessary to achieve an acceptalbe deconvolution result. All of this is accomplished in an iteractive mode of operation. First, the operator selects the offset correction, x(l), and the filter applied to x(k) and

y(k), and then judges the quality of the displayed deconvolution result. The step by step procedure is

as follows:

1. Preparation of the waveforms, x(k) and y(k).

a. Offset correction

b. Choose x(l)

2. Filtering

a. Choose the three filter parameters, L, S and A

b. Apply the filter to both waveforms, x(k) and y(k)

3. Deconvolution

a. Is the result stable, i.e., does it converge?

If "no" go back to step 1 or step 2.

b. Is the result smooth?

If "no" go back to step 2

c. Continue

4. Error Evaluation

a. Is the error acceptable?

If "no" go back to step 2

b. Continue

5. Print out results: d(k) with related quality parameters (mean error, standard deviation, etc.).

d(k) is the estimate of the impulse response.

The development of the appropriate filter (in 2a above) should nominally adhere to the following

rules

1. Start with the filter L = 2 , S = 3, A = 100

94
2. Keeping L = 2, and A = 10, Increase S until the result is stable (converges).

3. Keeping S and A constant increase L until the result is smooth.

4. After L and S are determined, compute the error criteria (section 3.4).

- if the mean value of the error, e, is greater than zero, decrease the value of A.

- if the mean value of e is less than zero, increase S until e becomes greater than zero. The smaller the value of e, the better is the result.

5. When e is close to zero, compute the standard deviation, 0. - if o is too large, increase L or S and adjust A to reduce O.

In summary, as it is evident from the above described procedures, time domain deconvolution is a progressive process. The programs used to implement time domain deconvolution were written to work in an iteractive mode. Knowing that the result must be convergent, smooth, and close to the true solution

(small error function), the operator selects the offset correction, x(l) , and the filter, judges the quality of the displayed result, and then modifies those selections until an acceptable result is obtained.

5.4.3.2 Estimation of the Insertion Impulse Response

of 300 Meters of RG 58 C/U Coaxial Cable

The following examples of time domain deconvolution employs the same x(k) and y(k) data as were used

in the frequency domain deconvolution experiments in sections 4.4.2 (two parameter method) and 4.6.3.2

(one parameter method). Consequently, the time domain results can be compared to the two frequency

domain methods, and vice-versa.

The input and output waveforms x(k) and y(k), respectively, are shown in figure 5.14; these are

the signals without additional noise added to them. Also, the time scales in figure 5.14 are expanded

(50 ns./div.) as compared to those in figures 4.9 and 4.10 (200 ns/div.).

Figure 5.15 shows the estimated impulse response d(k) using the filter parameters L = 2, S = 5, A = 3.15 applied to x(k) and y(k). Compare this result with those of figures 4.7b and 4.12 (2-parameter

F. D. method) and figure 4.34 (1-parameter F.D. method).

Figure 5.16 shows the estimated impulse response d(k) using the same filter L=2, S=5, A=3.15,

but with noise added to x(k) and y(k) to yield SNR-k and SNR-y values of 40 dB each. Compare this

result with figures 4.7f and 4.23 (2-parameter F.D. method) and figure 4.35 (1-parameter F.D. method).

95
h(k)

x(k)

i

I

Starts at k = 0, duration = N

^ H

»_ • I 12 3 n-2 N-1

Starts at k = 0, duration = M

• •

» I I i n* 12 3 M-2 M-1

No. of interior points

(points interior to zero- valued end points)

, = N-2

No. of interior points;

M = M-2

X

y(k) = x(k) * h(k)

Starts at k = 1, duration = L = N + M-1 Irstarcs at K = ±, a

I t f f-H—H J—

f

12 3 4 5

(a) Convolution

L-1

No. of interior

points

L = N + M-3

y

I I

h^(k) =

y(k)(l/'V)x(k)

Starts at k = 1, duration = N

I..

<—Hf-*-H^-«-^V4--^

12 3

(b) Deconvolution

N-1 N

No. of Interior

points: N, = N-2

Figure 5.1 The time positions and durations of x(k) , h(k) , y(k) and h^ (k)

.

x(k) and h(k) are two given sequences, each commencing at k = 0;

their durations are N and M, respectively. y(k) is calculated

by convolving x(k) and h(k), and the result commences at k = 1.

h^ (k) is h(k) calculated from y(k) and x(k) using deconvolution.

The calculations employ eqs. (5.2.1-8) and (5.2.2-6) or (5.2.2-21)

for convolution and deconvolution, respectively.

96
CO

I>Q- o>coo

<u

1-1 f^

1-1

C

<Uco

oo en \IO C g cn ^la 1-1 1 CO (U

iiJS 3p

*-« •H

x:H

X

CO rHd

00

<u

cfl •HOc

<u x;

4-1 Bo

v^

<4-l o

•u

0)e d

60

•H

97
o>cou

ree

60

COO

U3 ^ W ^—

rt ^>^

rH c

(1) to

4-1 <U /^ O ^—

to ^X

-a en c

•HO

ex:co

ei- m

(UVI

(U

tn&e

x;

toc3

Cu3O

-Oc

to3C

J-l D.3

V1-1 Ooo

4-1 e 3

00

1-1O

)-i toeou

toc

00

•H

CO x:

4-1 o

ja

)-i o

0)M3

00 H

98
o

>

o

u

0)

>

-ri U

u

0)

o

o

a:

T3 TJC

01 rt

>. ^

o

OJ

toc

o

a.

to

01

<i) to3

0)3

c

•H

tn

OJ

QJiJ

1a;to •Hoc

OJ

O

J-i OJe
Z>(f> o>coo

-o

0)>

(U

COo

03co

r4J

C3

0)

4-1

QJ

O

(-1 o o o

CO r-l

-a -ac

0) CO co c

0)I oo 0) Co oc

Qi: QJ ^UJ OJ £ 1 1- Vj AJ

en

a. Eou

CNin mcd E TJ 0) o^^ •H O

x;

)-4 -J »^ QJ u 00

100
co

•H

4-> 3

iHO>coo

0) -a

03C

paT3o

CO a:

i-j ^^c

0) CO

CJ ^-^

en ?X o en

0) 6ou

TDO

JZ

4-1 0)e I LO

OJM

00

•H

0)

J-i 3

oo

OO) 4-1 CO

U

"'g o

to

•H • ISIi-i •H cn cn CtKL

4-13 o

•H oo: s: Cu 4J xm 4J

-1 CO H- >^ O (U v-« ^ 4-) t3 •H

T3K

O a> CO

O

to •H

CO 4J3

4-1

a

J^ c • o T-l r-l

m CO Q y—S

t: c ^« ^ OJ

4-1 •H i-iO ex: a, m \

0)

(n 3

"a ^1- o oo

CO ^. en QJ ^1 2 4-1 o

101
01N

•Hu

a;

o

M

Co

cn0)

3

2:yj<r

en

OJ

U

3

00

noO

x:

4J

01e

c

o

o

>

c

o

<J

102
T3

(UN

•HU

CO tH3

00

0)u

JO

•acO

COco

a.

en

tu

^1

q;

en iH3a

4)3

GOO

4J

0)

cao<eooo

4-1 0)eco

'^ o ^ >

13

OO

0) 0)

J-1 -a

COe oj

•U ^

tnu

0) nj 3

00 •H

103
I

I

(a) L- Filter. L = 3,

f = 340 MHz. Vertical

CO

scale 10 dB/Div. Horizontal

scale 50 MHz/Div.

CO

Pi(b) S, A-Filter, = 5, A = 1

f = 250 MHz. Vertical

CO

scale 10 dB/Div. Horizontal

scale 50 MHz.

CO

(c) Combination of L and S Filters. The L filter

f = 340 coincides with

CO

the S, A Filter first maximum.

Vertical and horizontal scales

are 10 dB/Div. and 50 MHz, Div.

Figure 5.9 Transfer properties L and S,A filters and their cascade properties,

104
A = 0.50

A = 1.00

A = 1.50

A = 1.80

A = 2.00

A = 2.05

A = 2.08

A = 2.40

A = 5.00

Figure 5.10 Change in the transfer properties of the S, A

Filter as A is increased. S = 5.

105
^ ^

L5 59

L5 59

A3. 85

\

1 A3. 85

\\

\ /^ ^ \J\ \ / -N

\(] \ /

i

\

, ,/

\

L5 58 A2.95

L5 58 A2.95

\,

N

\

\V /

\ \^~N /

1 r

> /

/ li \

1 "X L4 48 A3. 35

\,

L4 58 A3. 35

k

\

\ <^~- \

\ r^ nl^ \

f \ i

' \

J

i

^'s. L5 57 A2.30

\

L5 57 A2.30

\

\

\ A

\ / \

^

i f V, (

"^s

1 \

) \

1

1 (

"^^ L4 57 A2.95

>N. L4 57 A2.95

\

\

\ / ^^

\ /i \ r \

\\

-\ L3 56 A2.40

\,

L3 56 A2.40

\

\

\ /I

rxkA r\^-^

1 V \

\

j

1

--^ L3 55 A1.50

X L3 55 ^1,50

s,\

\

/\

\

/ \

1

,

^

1

Vertical scale 10 dB/Div. Horizontal scale

Af/Div. is a function of the time window T, Af = 1/T.

Figure 5.11 Examples of composite filters.

106
"N iv L 8 51^ ^ 5- \

\

\

' A/\ l\

\

1 \

\( rr\ o-

^

\i

V

|-

1 \

1 i

S

N J -1 Ri/, M^ on \

\

\

r\r^^\'\

\ ^ ,/\ /

6

1^ \l hir^

i i 11

1

^ N T,7 SI 3 A 5.60

\

^

\

\

«

\ /\/

^

V

\ \ \

/ K r

1 \ \

>

1 ll \/

1

^~^ N,

\

L6 512 A5.15

\

\

\

/[\r \ r^

1

^-> r^ N

\l

-\ L5 511 A4.70 N L 5 '510 AA.in

\

V

\

\

/

\j ' S

\ — -vl=» -—/ =v- {\f\ /^ y A ^

1 \ \

\J \ ]

^^-^ r J

I \

1

\

]

1

—-- N

\, L6 59 A3. 25

\

\

I M \ ^ -s.

'\ \{

•^ ^ H[/ ^

\\ \ 1 ^

1

1 \

1

I

1 1

N L5 59 A3. 80

\

\

^

\ /^ I Ni/'^ \ -\

V 1 \

1

1

'

Figure 5.12 Additional examples of composite filtering scales are

the sam.e as in figure 5.12.

107
O 3 CO (U 60

^ CD c c o c

ro

1

CO S-i o -H • J-) -H cn <

Q O •H g .-V j_i ;.j .

11

Vj -u 4J M CO -T3 o *H <r O O (U ^ OJ 0) CO M 14-1 • 00 (U 4-1 (1) CO ^ (X

1 o c tn <u 00 3 a; cn

en o -H M -a 0) CO QJ .H • i (U iJ-i 13 OJ 4-J C -H >^

t d c 4JtnC(U "^ojiH M .-1 - o COJ-HXfU'HS-ia) t ca >^ cu •H d^— j-ii—1 en co >

1

1 > o en •H CL tn -H C OJ a)^cncu6o<4-i < 0) ^J jz 0) v-i CO a. cj "

1

r-\ 3 i-i H iH 3 CO X -73 (U ITI -n cr o CO -H 0) O C Cl< 1

• C QJ o J-I CO cn S CO CO

\ '

CN n3 !-i d > <u s-1 -u <u \

U-l 0) O O > O <Xi u

9 C/3 C a N c/2 M fiH i

U-i -H m -H 3 N Sj *v \ U-l i« iH ^J 0) a B o .H O \ O O T3 O iH • S '« ^

Ln ^ »

1 -H (U ^ ^ T3 < • Csl C J-J CO Vj -H QJ O CO LO O 3 O •H (U tn X OJ <N OJ •H CJ -H tn X tn CO j-i .H 3 *s U XJ QJ W O CO CO iH m

m O C )-i -3 O. ID 6 II CO <H

\ ' 0) 0) CU JZ -H >

iTl .H > > QJ 4-1 QJ OJ iJ O T3

V I vO OJ -H jd-HX^tn GM-iC

\ 0^ CO M CO 4-IS4-1HQJM.IOC0

\ • \ < n

\

\

> CO

\\ oo"

\ LO

y-i \ \

(U

^

• V \\

3

00

>. i •H u

\

fe C \ o . o

^^

\ ^"

cr

0)

J-i

CN

1

< 0-) h " <

\ •

M-l \

4-( \

\

O

1 "

^ 00 •

J-1 3 \

\

\

< irn O

\

v O

"

. \^ CN \1

1

( »^" "

o ^" ^" i_ ^"^^ *r

r-l \ 1

I s \ ro - \

\ \

\ r-. ^

«

- N

\

\

\

r-i O

CN " .... * -^ •* * * ^ ^ " V

L 1.

^1

Lno

C7^•

O _ ^^^! -. ^v N

> ^O O *^ "^^^

i %

c^ ~3- LO vD

II II II II < <: < <

V3 LTi vO r^ c30 cn o H ^ 0-1 ^t Ln

iH rH 1—1 .H t-H rH

108
IhFUT UAUEPQRH

Diu e.TeeE i mx « 8.638E 2

HiN —e.Teee i

HQRIZONTM. SCALE

DIU = SB.ee NS

SIO«L TO NOISE

75.1 OB

OUTPUT WWEFORM

DIU = e.seeE e

MAX « 8.45eE I

MiN —e.seee e

H0RI20KTM. SCfiLE DIU « 9e.ee hs

SICH^ TO NOISE

46.5 06

t\

\

\\

\

\V"^

Figure 5.14 Input and output waveforms x(k) and y(k) used in the frequency

and time domain deconvolution experiments before noise is added

to the signals. These waveforms are identical to those in

figures 4.9 and 4.10 except that the abscissae are 50 ns/Div.

(as compared to 200 ns/Div.).

109
RESULT CF TX;<E DOMAIN OECONUXUTION

OIU » Q.250E. -Z

f%5X « e.22SE -1

WIN —e.25ee -2

HORIZONTAL SOALE

DIU > 58.88 HS

FILTER « L 2.S 5. A 3.13

ERROR FU^TION

f\

\

I

\

\j

V

\\ V -

OIU = 0.547E -2

HfiX = 8.26SE -1

MIN "-8.278E. -1

.Tifllll , 1 liii .1 .1 J . 1, 1 . .1 f€AM = e.214E -6

SI(»^ 8.eB6E -2

HORIZQKTAL SCdJE.

OIU - 98.88 NS

FILTER ' L 2.S 5. A 3. IS

r ilel(FJ!Mlir'i^!!f!I!l^^F.m'

|[ill: MIOIi'iLiiUiiliJiiiyililiiiblki

1 '1 11 1 1

i

1 1

1 II

Figure 5.15 Result of the time domain deconvolution d(k) the L2 and S=5 , A=3.15

filters applied to the same waveforms used in the frequency domain

deconvolution experiments (figure 5.14). Compare to the frequency

domain deconvolution result, figure 4.12.

110
l^eSULT OF Tire OOmiN OECONUOLUriON

DIV « e.258E -2

i\

flAX « 8 22SE -1

MIH —e.25BE -2

\

\

\

K

HORIZONTfiL SCfiLE \1 OIU « 9B.ee KS

V Kynh^ \A« 4

FILTER ' L 2.S 5. A 3.ie

m.KA

-"f^ m{ -U^ M^^M ^

E?^=:::? fjnction

oiy « e.i63E -I WW = e.8l8E -1

WIN «-e.835E -1

rtAN « e.676E -4

siGw«= e.3ieE -1 nimTi"^^!

HORIZOHTflL SOLE J. RJI'iPlf ffH^if fll

Diu - 5e.ee hs ::: ::_:'^":^t|i

L 2.S 5. ft 3.18 j

Figure 5.16 Result of time domain deconvolution d(k) when the input and output

waveforms are corrupted by the addition of 35 and 6 dB of noise,

respectively, (figures 4.20 and 4.21) to yield SNR-x 40 dB and

SNR-y 40 dB. The filters applied to each waveform are the same as those in figure 5.15. Compare this result to the frequency domain

result, figure 4.23.

ill
6. SUMJIARY AND CONCLUSIONS

6 . 1 Sununary

The objective of the studies reported here was to investigate methods of deconvolution to be applied

to time domain sequences such as those that represent the signals in sampled or discrete data systems.

In the preceeding chapters several deconvolution methods have been explored analytically and experimentally. After some introductory comments and some general presentations of the deconvolution problem

and error criteria, two distinct domains for deconvolution of time domain sequences were investigated:

(1) the frequency transform sequence domain, or frequency domain, and (2) the time sequence domain, or time domain.

Two methods of discrete frequency domain deconvolution were studied; the first used a two parameter

filter whose parameters were subjectively selected by the operator to filter the transfer function H(n)

prior to its inversion to the time domain to yield the estimated impulse response, d(k). The second

method was one-parameter optimal technique in which the frequency domain filter was optimally adjusted to minimize an objective function consisting of the weighted sum of an error-energy and roughness-energy

(departure from smoothness). The two methods were applied to the measurement of the impulse response of

a 300 meter long coaxial cable for which experimental input and output time domain sequences were known; various degrees of computer generated noise (pseudo-random sequences) were added to the relatively noise

free experimental data sequences in order to simulate the deconvolution of noisy data. Also, in addition

to the coaxial transmission line experiments, the two-parameter method was applied to the measurement of the impulse response of a broadband antenna, while the one-parameter method was employed to obtain the impulse response of a waveshaplng filter.

In the time domain studies, two methods were also investigated, discrete classical deconvolution

and discrete iterative deconvolution; however, the latter method was not applied to any experimental data. The classical method directly inverted the discrete convolution equation to obtain the unknown impulse

response. In the iterative method, the discrete convolution equation is used to form successive approximations to the impulse response; the resultant approximate impulse responses are convolved with the input signal to yield an output sequence which is then compared to the known output sequence. When the

two sequences agree within some specified limits, the iteration process is stopped and the N-th iterated

impulse response is declared to be equal to the sought-after or unknown impulse response. Both methods were applied in a simulation study using specified x(k), h(k), and y(k) sequences uniquely related to each other through the convolution equation, in the absence and presence of simulated noise. Also, a simple regularization filter was used to demonstrate the improvement provided by filtering to offset the effects of noise. The classical method and a more sophisticated (composite) filter were applied to the same coaxial cable experimental data which were used in the frequency domain deconvolution studies;

thus the results of the frequency and time domain methods could be compared to each other. The synthesis and properties of the composite discrete time domain were discussed and developed.

6.2 Conclusions

For the deconvolution methods considered here the frequency domain one-parameter method provided

the smoothest results without any obvious aberations such as Gibbs-ef f ect. The two-parameter method

may be useful in some interactive filtering situations, but for the most part, it is only of academic

interest. The one-parameter method is an optimal filtering method, and as such the weighting between

error and smoothness (energies) can be controlled.

The time domain deconvolution methods are still in their infancy and are purely experimental.

Consequently, much work remains to be done to bring them to fruition so that they can be easily applied.

112
Nothing encountered in this investigation rules out the promise held for time domain methods, i.e.

the elimination of transform filtering.

113
7. REFERENCES

[1] McGillem, C. D.; Cooper, G. R. Continuous and discrete signal and system analysis. New York:

Holt, Reinhart and Winston; 178; 188; 1974.

[9

[10

[11

[12

[13

[14

[15 [16

[17 [18

[19

[20

[21

[22 [23

[24

[25

Brigham, E. 0. The fast fourier transform. New Jersey: Prentice-Hall; 1974. Tikhonov, A. V.; Arsenin, V. Y. Solutions of ill-posed problems. Washington, D.C.: 1977. Biraud, Y. G. Les methodes de deconvolution et leurs limitations f ondamentales . Revue de physique

appliquee. Tome 11, No. 2: 203-214; 1976 March.

Jones, A. F. ; Misell, D. L. The problem of error in deconvolution. J. Phys. A., Gen. Phys. 3: 462-472; 1970. Hildebrand, F. B. Introduction to numerical analysis. New York: McGraw-Hill; 1956. Von Blaricum, M. L. ; Mittra, R. Problems and solutions associated with Prony's method for processing transient data. IEEE Trans. EM Compat. EMC-20(1): 174-182; 1978 Feb. Riad , S. M. ; Nahman, N. S. Application of the homomorphic deconvolution for the separation of TDR signals occurring in overlapping time windows. IEEE Trans. Instrum. & Meas. IM-25(4): 388-391;

1976 Dec. Cans, W. L. Present capabilities of the NBS automatic pulse measurement system. IEEE Trans. Instrum. Meas. IM-25(4): 384-388; 1976 Dec. Cans, W. L. ; Andrews, J. R. Time domain automatic network analyzer for the measurement of RF and microwave components. Nat. Bur. Stand. (U.S.) Tech Note 672; 1975 September.

Op. Cit. [1]. Section 5-18.

Op. Cit. [2] . Chapter 8. Digital Signal Processing Committee, ed. Programs for digital signal processing. New York: IEEE

Press; 1979. Oppenheim, A. V.; Schafer, R. W. Digital signal processing. New Jersey: Prentice-Hall; 1975.

Op. Cit. [1]. Section 5-15.

Op. Cit. [13]. 405, 410. Kuo, Franklin F. Network analysis and synthesis. New York: John Wiley and Sons; 1962; Chap. 13. Twomey, S. The application of numerical filtering to the solution of integral equations

encountered in indirect sensing measurements. J. Franklin Inst. 279: 95-108, 1965 Feb. Hunt, B. R. Deconvolution of linear systems by constrained regression and its relationship to the Weiner Theory. IEEE Trans, on Automatic Control AC-17(5): 703-705: 1972 October.

Op. Cit. [14]. 204. Personal communication from Dr. Yves Balcou, Dept. of Crystalline Physics and Chemical Structure,

University of Rennes, Rennes , 35000, France.

Op. Cit. [14], p. 125.

Op. Cit. [1], p. 37 Cans, W. L.; Nahman, N. S. A note on the conversion of an experimentally obtained step-like waveform into a duration limited waveform for implementing Fourier transform analysis. Submitted for publication IEEE I&M Transactions.

Andrews, J. R. ; Pulse reference waveform standards development at NBS. Proceedings of the ATE

Seminar/Exhibit, Pasadena, CA. , Jan. 1981, pp. IV-13-IV-19. Benwill Publishing Corp., Boston, MA.

02215.

Op. Cit. [1], p. 89. Janson, Peter A. Method for determining the response function of a high-resolution infrared

spectrometer. J. of the Opt. Soc. Amer., Volume 60, No. 2., 194-191, Feb. 1970.

114
[28] van Cittert, P. H. Z. Physik, Vol. 69, 298, (1931).

[29] Op. Cit. [1] , Chapter 7.

[30] Pipes, L. A.; Harwill, L. R. Applied mathematics for engineers and physicists", 3rd Ed., pp. 562- 567, McGraw-Hill Book Co. 19 70.

115
NBS-n4A iREv. 2-8C) U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO. NBS TN-1047

2. Performing Organ. Report No 3. Publ ication Date October 1981

4. TITLE AND SUBTITLE

DECONVOLUTION OF TIME DOMAIN WAVEFORMS IN THE PRESENCE OF NOISE

5. AUTHOR(S) Norris S. Nahman and M. E. Guillaume

6. PERFORMING ORGANIZATION {If joint or other thon NBS. see ;n struct/onsj

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

8. Type of Report & Period Covered

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State, ZIP)

10. SUPPLEMENTARY NOTES

[2J Document describes a computer program; SF-185. FlPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most si gnificant information. If document includes a significant

bibliography or literature survey, mention it here) Deconvolution or inverse filtering was used to determine the impulse response of

a system using noisy input and output time domain sequences (discrete data) . Frequency

and time domain methods were studied along with the synthesis of the filters required

to obtain stable and smooth results. For the methods studied it was concluded that the

superior technique was provided by an optimal frequency domain method implemented via

the FFT. Also, it is pointed out that the time domain methods are only in their

infancy and still retain the promise of avoiding transform domain filtering. Examples

are presented in which the impulse responses are determined in the presence of varying

degrees of noise for a coaxial transmission line, a wave-shaping filter, and a broadband antenna.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)

deconvolution; impulse response; inverse filtering; noise contamination ; time domain;

waveforms.

13. AVAILABILITY ^^ Unl imited Q For Official Distribution. Do Not Release to NTIS

[231 Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

20402.

[X] Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

122

15. Price $11.00

USCOMM-DC 6043-PPn

<IU.S. GOVERNMENT PRINTING OFFICE:1981—576-001 / 1206
NBS TECHNICAL PUBLICATIONS

PERIODICALS

JOURNAL OF RESEARCH—The Journal of Research of the National Bureau of Standards reports NBS research and develop- ment in those disciplines of the physical and engineering sciences in which the Bureau is active. These include physics, chemistry,

engineering, mathematics, and computer sciences. Papers cover a broad range of subjects, with major emphasis on measurement

methodology and the basic technology underlying standardization. Also included from time to time are survey articles on topics

closely related to the Bureau's technical and scientific programs. As a special service to subscribers each issue contains complete

citations to all recent Bureau publications in both NBS and non- NBS media. Issued six times a year. Annual subscription: domestic

$\}: foreign $16.2.'^. Single copy, $} domestic; $3.75 foreign. NOTE: The Journal was formerly published in two sections: Section A "Physics and Chemistry" and Section B "Mathematical

Sciences." DIMENSIONS/NBS—This monthly magazine is published to in- form scientists, engineers, business and industry leaders, teachers, students, and consumers of the latest advances in science and

technology, with primary emphasis on work at NBS. The magazine

highlights and reviews such issues as energy research, fire protection, building technology, metric conversion, pollution abatement,

health and safety, and consumer product performance. In addition, it reports the results of Bureau programs in measurement

standards and techniques, properties of matter and materials, engineering standards and services, instrumentation, and

automatic data processing. Annual subscription: domestic $11;

foreign $13.75. NONPERIODICALS

Monographs—Major contributions to the technical literature on

various subjects related to the Bureau's scientific and technical activities. Handbooks—Recommended codes of engineering and industrial practice (including safety codes) developed in cooperation with in- terested industries, professional organizations, and regulatory

bodies. Special Publications—Include proceedings of conferences sponsored by NBS, NBS annual reports, and other special publications

appropriate to this grouping such as wall charts, pocket cards, and

bibliographies. Applied Mathematics Series—Mathematical tables, manuals, and

studies of special interest to physicists, engineers, chemists,

biologists, mathematicians, computer programmers, and others engaged in scientific and technical work.

National Standard Reference Data Series—Provides quantitative data on the physical and chemical properties of materials, compiled from the world's literature and critically evaluated. Developed under a worldwide program coordinated by NBS under

the authority of the National Standard Data Act (Public Law

90-396).

NOTE: The principal publication outlet for the foregoing data is the Journal of Physical and Chemical Reference Data (JPCRD)

published quarterly for NBS by the American Chemical Society

(ACS) and the American Institute of Physics (AlP). Subscriptions,

reprints, and supplements available from ACS, 1 155 Sixteenth St., NW, Washington, DC 20056.

Building Science Series—Disseminates technical information developed at the Bureau on building materials, components,

systems, and whole structures. The series presents research results,

test methods, and performance criteria related to the structural and

environmental functions and the durability and safety characteristics of building elements and systems.

Technical Notes—Studies or reports which are complete in themselves but restrictive in their treatment of a subject. Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject area. Often serve as a vehicle for final reports of work performed at NBS under the sponsorship of other government agencies.

Voluntary Product Standards—Developed under procedures

published by the Department of Commerce in Part 10, Title 15, of

the Code of Federal Regulations. The standards establish nationally recognized requirements for products, and provide all concerned interests with a basis for common understanding of the

characteristics of the products. NBS administers this program as a supplement to the activities of the private sector standardizing

organizations. Consumer Information Series—Practical information, based on NBS research and experience, covering areas of interest to the con- sumer. Easily understandable language and illustrations provide

useful background knowledge for shopping in today's tech- nological marketplace. Order the above NBS publications from: Superintendent of Documents, Government Printing Office, Washington, DC 20402. Order the following NBS publications—FIPS and NBSIR's—from

the National Technical Information Services, Springfield, VA 22161

Federal Information Processing Standards Publications (FIPS

PUB)—Publications in this series collectively constitute the Federal information Processing Standards Register. The Register

serves as the official source of information in the Federal Govern- ment regarding standards issued by NBS pursuant to the Federal Property and Administrative Services Act of 1949 as amended.

Public Law 89-306 (79 Stat. 1127), and as implemented by Executive Order 11717(38 FR 12315, dated May II, 1973) and Part 6 of Title 15 CFR (Code of Federal Regulations).

NBS Interagency Reports (NBSIR)—A special series of interim or

final reports on work performed by NBS for outside sponsors

(both government and non-government). In general, initial dis- tribution is handled by the sponsor; public distribution is by the National Technical Information Services, Springfield, VA 22161,

in paper copy or microfiche form.
U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Washington. D.C. 20234

OFFICIAL BUSINESS

Penalty for Private Use, S300

POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCE

COM-215

SPECIAL FOURTH-CLASS RATE

BOOK