&EPA
United States
Environmental Protection
Agency
Office of Radiation Programs
Radiation
Eastern Environmental
Radiation Facility
P 0 Box 3009
Montgomery AL 36193
EPA-520/5-80-003
March, 1980
Radiological Data Analysis
in the Time and
Frequency Domain
Prepared for EPA by
AUBURN UNIVERSITY AT
MONTGOMERY UNDER EPA
CONTRACT No. 68-02-3049
-------
RADIOLOGICAL DATA ANALYSIS IN THE
TIME AND FREQUENCY DOMAIN
by
Dr. Donald A. Chambless
Department of Mathematics
Auburn University at Montgomery
Montgomery, AL 36117
Contract No. 68-02-3049
EPA Project Officer: Dr. Jon A. Broadway
Prepared for
Environmental Protection Agency
-------
DISCLAIMER
This report was furnished to the Environmental Protection Agency by
Auburn University at Montgomery. Montgomery, Alabama, in fulfillment of
Contract No. 68-02-3049.
The contents of this report are reproduced here.
in as received from the contractor.
The opinions, findings, and conclu-
sions expressed are those of the author and not necessarily those of the
Environmental Protection Agency.
Mention of company or product names is
not to be considered as an endorsement by the Environmental Protection
Agency.
-------
TABLE OF CONTENTS
Page
I. I NTRODUCT I ON. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
I I. SUMMARY.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4
I I I. CONCLUS IONS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .11
IV. RECOMMENDATIONS..............................................13
V. SECTION 1: NOISE REDUCTION VIA FREQUENCY
DOMAIN ANALYSIS................................. .16
VI. SECTION 2: GENERAL CONSIDERATIONS IN
RESOLUTION IMPROVEMENT...........................29
VII. SECTION 3: QUADRATURE METHODS AND SINGULAR
VALUE ANALYSIS.................................. .32
VIII. SECTION 4: GENERAL PHILOSOPHY OF REGULARIZATION
METHODS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .35
IX. SECTION 5: DISCRETE IMPLEMENTATIONS OF
REGULARIZATION. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .37
X. SECTION 6: CONTINUOUS/DISCRETE REGULARIZATION...............43
XI. SECTION 7: THE PROBLEM OF D.L. PHILLIPS.....................49
XII. SECTION 8: APPLICATION TO Ge(Li)
SPECTROMETER DATA................................ 92
XIII. BIBLIOGRAPHY................................................ .138
-------
INTRODUCTION
The accurate measurement of levels of environmental pollutants is a
matter of central importance to the United States Environmental Protection
Agency.
However, the level of radioactivity encountered in samples of
man-related nuclide pollutants is typically exceedingly low, and the highly
sensitive state-of-the-art instrumentation employed in detection worK be-
comes (because of its very sensitivity) affected not only by the radiation
of interest but also by many other background radiation effects, as well
as random effects.
Currently, then, it may often be desirable to apply
data enhancement techniques (of an analytical nature) to the raw instru-
mentation data.
There are two types of enhancement of the data which may be considered
in any given application.
First of all, any instrument data record typi-
cally shows the presence of two signal components:
a deterministic com-
ponent, which reflects the "genuine" informational content of the data
record, and a random component which results from stocha~~ic behavior of
the basic phenomenae.
Usua lly one woul d 1 i ke to separate these" signa 1"
and "noise" components and retain only the former.
It is often possible to
partially achieve this objective using methods described in Section 1.
Secondly. even assuming that a noise-free data record from the analyzer
has been produced, there still may be an unsatisfactory level of resolution
of the recorded signal due to simple inadequacy of the resolving capacity
of the detector and instrument system.
In such a case, if quantitative
measurements are to be obtained then resolution improvement of the record
through analytical means must be achieved; this is, both in principle and
practice, a much more difficult proposition, as the details of this report
1
-------
will amply demonstrate.
This report contains an outline of work accomplished under EPA
Contract Number 68-02-3049 during the year beginning with September 27~
1978.
Under the provisions of this contract~ into which Auburn University
at Montgomery and the Eastern Environmental Radiation Facility ("EERF") of
the U. S. Environmental Protection Agency entered~ the author was con-
tracted to devote approximately two-thirds of the year mentioned to an
investigation of the characteristics of gamma radiation data from Ge(Li)
spectrometers and the applicability of data enhancement techniques (of
the two types outlined above) to this data.
Under the provisions of this
contract the author travelled the short distance from the university to
EERF after completing classes~ office hours~ etc.
Thus the actual contract
work was largely accomplished at the laboratory site~ with immediate access
to not only the instrumentation of interest but also the computing and
library facilities of EERF and~ most importantly~ the continual coordination
and supervision of the Project Officer.
It is felt that this was a most
efficient manner of administering the project~ both in terms of minimizing
costs and maximizing the effectiveness of the technical supervision.
In
particular~ the contract overhead rate was only twenty per cent.
Work done jointly by the author and the Project Officer during their
investigations related to the contract research topic have resulted in a
number of accomplishments.
These are summarized in the following:
(1)
D.A. Chambless and J.A. Broadway~ Fourier Transform Methods for
Analysis of Ge(Li) Spectral Data~ Presented to the Health Physics
Society~ Atlanta~ GA~ July 6~ 1979.
(2 )
D.A. Chambless and J.A. Broadway~ Digital Filtering of Speckle
2
-------
Interferometric Data, Presented to the Society for Experimental
Stress Analysis, San Francisco, CA, May 23, 1979.
(3)
Jon A. Broadway and Don A. Chambless, Spectral Resolution Improve-
ment via Regularization Methods and Singular Value Decomposition,
Presented to the Health Physics Society, Philadelphia, PA, July
11, 1979.
(4)
D.A. Chambless and J.A. Broadway, Digital Filtering of Speckle-
Photography Data, Published in Experimental Mechanics, Vol. 19,
pp. 286-289, August, 1979.
(5)
J.A. Broadway and D.A. Chambless, Resolution Improvement of
Gamma Radiation Spectrometer Data, Presented to the International
Symposium on Ill-Posed Problems:
Theory and Practice, Newark,
DE, October 2, 1979.
(6)
J.A. Broadway and D.A. Chambless, Constrained Regularization
Methods for Spectral Resolution Improvement (manuscript in
preparation).
Finally, it should perhaps be noted that the original proposal sub-
mitted contained also a proposal for an optional one-year extension of the
contract effort.
At this writing, the Office of Radiation Programs has
exercised its option and effected such an extension.
Therefore, the work
begun under this contract arrangement will be continued for another two-
thirds man-years and a report pursuant to that work and extending the work
reported herein will be forthcoming in September 1980.
3
-------
SUMMARY
As mentioned in the Introduction, the work conducted in this contract
effort consisted primarily of feasibility studies concerned with the reduc-
tion of noise levels and the improvement of resolution of radiation spectra
recorded by Ge(Li) gamma spectrometers.
The methods investigated in the
course of this work included a number of techniques already recognized in
the literature as well as some not yet appearing (to the knowledge of the
author) .
Closely related problems in other fields ranging from photo-
graphic enhancement to forensic science to historical document restoration
have been considered in certain instances.
These efforts were made so as
to allow comparison of the methods considered under this contract activity
with other existing methods and to provide a background of results with
"well-known" problems for consideration.
The noise reduction problem is considered in Section 1 of the main
body of this report.
In this section the foundations of a frequency domain
approach to noise suppression are described.
Under the assumptions con-
cerning the data detailed in this section, considerable increase of signal-
to-noise ratio is often achieved through the application of the methods
described, especially when the deterministic component of the composite
signal represents a function which can be well approximated by low order
trigonometric polynomials.
On the other hand, if a signal should be
qualitatively characterized as, say, the superposition of a linear back-
ground and a low order trigonometric polynomial plus noise then one would
expect that the removal of the linear trend (through standard least squares
methods) would be necessary as a preprocessing procedure before the filter-
ing methods of Section 1 are applied; this expectation is based on the fact
4
-------
that the linear carrier signal is not at all well approximable through
lower order Fourier expansions.
In any case, the methods of Section 1 are found to be very inexpensive
to implement due to the availability of very efficient ("fast") Fourier
transform algorithms.
In actual application, the speed and convenience of
the use of the techniques are very dependent upon access to efficient
graphics equipment and software.
Unhappily. there has been no opportunity
for such use of high-speed graphics during the course of this study, and
so the number of applications actually included herein is rather small.
Nonetheless, some rather striking results were obtained, as a comparison of
Figures 1.1 and 1.5, for example, reveals, and the quality of these en-
hancements has been recognized through the publication of [6].
In Section 2 the mathematical generalities concerning the problem of
attempting to analytically enhance an analyzer data record by correcting
for limitations in instrument resolution characteristics are given and the
specific conceptual problems of interest are stated.
In Section 3 the
futility of straightforward quadrature approaches to this exceedingly
difficult problem is suggested (although specific indication of this is de-
layed until Section 7), and the singular value decomposition of a matrix
is formulated and its importance discussed.
In particular, the singular
value analysis approach to the solution of ill-posed linear problems,
which has found considerable popularity among those investigators consid-
ering the processing of photographs [23,24], is motivated and described.
The major concentration of effort in this contract endeavor came to
be the study and application of the family of "regularization" methods
[22, 27-30] for the solution of ill-posed linear operator equations.
The
5
-------
general conceptual framework of regularization methods is given in Section
4, and a number of discrete implementations of regularization are formu-
lated in Section 5.
All of the methods of Section 5 are basically matrix-
theoretic.
Some employ the difference operator in order to allow (discrete)
considerations of the first (or higher) derivatives of the functions of
interest to be included while the most basic form of discrete regularization
does not.
Later (in Section 7) the effect of this higher generality is
examined.
Although there appears to be some advantage in including consid-
eration of one or more derivatives in regularization calculations, there
was not a generally marked improvement in the quality of the results noticed
while, on the other hand, the cost of doing "differential regularization"
was significantly higher.
It is known that more efficient algorithms for
these advanced methods can be developed (with an additional expenditure of
time) but there seemed to be no justification for embarking on such an
effort at this point.
In Section 6 a form of regularization [32] ("CDR regularization")
which represents a very significant conceptual advance over the discrete
methods of Section 5 is discussed at some length.
During the course of
this contract year this was actually the first form of regularization
which was considered, and considerable time and effort were devoted to
this aspect of the project.
In particular, the method was discretized (as
any computational method must eventually be) and applied to model problems
such as those in Sections 7 and 8.
Results from these efforts showed a
striking improvement over the singular value analysis (and harmonic
analysis) solution attempts which had been previously conducted.
Similarly,
the subsequent discrete regularization calculations (using Simpson IS rule
6
-------
to perform quadratures as previously) were not nearly of the quality of
the results using the methods of Section 6.
Sometime later it was dis-
covered that Simpson's rule itself was causing significant degradation of
the results obtained through singular value analysis or discrete regular-
ization (as compared to those achieved by means of the simple rectangle
rule, for example). Of course, one ordinarily thinks of Simpson's rule
as being generally more accurate than the rectangle rule for most usual
types of functions, and hence an adequate explanation of this phenomenon
is not immediately apparent.
(An indication of the reason for this is
given in Section 7.)
In any case, the COR method was found to give results
of about the same quality with either quadrature rule and therefore appears
much more powerful (for the types of problems being considered) than the
methods of Sections 3 or 5.
A fairly extensive effort was made to lower
the initially high cost of the COR computations and a good deal of success
was achieved in this regard; details are given in Section 6.
Sections 7 and 8 are devoted to outlining an indication of typical re-
sults obtained with resolution enhancement problems during the course of
this contract effort.
Two basic types of problems were considered.
The
first problems, which are discussed in Section 7, are derivatives of the
original problem of D. Phillips which was considered in his 1962 paper [22].
In this model problem, which has come to be a classic example in the 1it-
erature of the field, the resolution improvement of noisy pulse-like data
is considered.
Since the data of this problem is somewhat typical of the
gamma spectral data of ultimate interest in this contract effort and since,
also, this problem of Phillips is one of the few prototypes of the field,
there was a considerable amount of effort devoted to "sorting out" the
7
-------
methodology considered during the contract year by means of analytical
investigation of Phillips' problem.
It was in the course of this work that
the sensitivity of the calculations to the selection of the quadrature rule
was noticed.
Clear indication of this phenomenon is included in Section 7.
It is of considerable interest to note that the degradation generally intro-
duced by Simpson's rule was largely overcome by the DER and COR methods
described in Sections 5 and 6, respectively.
On the other hand, all of
the methods considered (with the notable exception of the simple quadrature
and Fourier deconvolution methods) were found to be reasonably capable of
handling Phillips' problem when the rectangle rule is used for quadratures.
In Section 8 the second basic type of resolution improvement problem is
considered.
Here Ge(Li) spectral data from radium sources as recorded by
a Nuclear Data NDIOO spectrometer are investigated by means of the method-
ology of Sections 5 and 6 in two basic situations.
In the first group of
calculations the data used were realized by the superposition of 5% random
noise on the impulse response of the spectrometer (as defined by the data
record resulting from a monochromatic radium energy line); this was refer-
red to as the "spectral data". The other g ro u p of calculations was per-
formed using data from a different radium peak; this data was called the
"cross-spectral data".
It should be pointed out that although Ge(Li)
spectrometers are not generally thought of as shift-invariant devices, it
was felt that this was a reasonable assumption herein due to the small
energy ranges considered in Section 8.
vJith both of these sets of data it was found to be exceedingly diffi-
cult to produce a result sharply approximating the ideal discrete pulse for
which one might hope under perfect circumstances.
The instability of the
8
-------
basic problem per se is such that even a few per cent noise superposed on
the impulse response causes extreme ~c;"tJLAtational difficulty and results
closely approximating the ideal one seem impos~c~le to obtain without some
further modification of the problem statement beyond that provided by the
various regularization tec~niques.
Two such modifications were considered,
the first representing a
compromise" problem statement and the other
involving the imposition of constraints on the computed solution.
The compromise used was effected by selecting a Gaussian pulse with
half-width significantly smaller than that of the impulse response of the
analyzer and then modifying the problem statement such that the ideal solu-
tion of the modified problem* consists of this selected Gaussian.
Philo-
sophically, one may view this as solving for a certain moment of the solu-
tion instead of the solution, per se.
Thi s II data spreadi ng" method was
found to increase the stability of the calculations in Section 8 while not
causing any large increase in the computational costs.
The relevant corn-
parison can be drawn from Figures 8.8 and 8.11, for example.
On the other
hand, a second modification of the problem was constructed by explicitly
incorporating the a priori knowledge of nonnegativity of the solution.
Very striking improvements in the stabil ity and sharpness of the results
obtained were found to be gained through this type of consideration as
Figures 8.38 and 8.40 reveal.
There was no remaining time in which to
*C. Nelson of the Environmental Analysis Division, Office of Radiation
Programs, suggested that such a modification be considered.
9
-------
conduct further calculations of this type and thus spectra larger and m~re
complex than those indicated in Section 8 have not yet been considered at
this writing.
The conduct of such larger scale investigations will be an
early topic in the work of the year of continuation of this contract effort
10
-------
CONCLUSIONS
Two types of enhancement procedures for Ge(Li) gamma spectral data have
been considered under this contract effort:
noise reduction techniques and
resolution improvement methods.
In the case of noise reduction the results
of Section 1 (as observed, for example, by comparison of Figures 1.6 and 1.7)
show that very significant increases in signal-to-noise ratio can be obtain-
ed by means of fast (and inexpensive) frequency domain procedures (under
certain hypotheses concerning data and noise characteristics).
One does
have to make a significant effort in software development in this regard
and high-speed graphics equipment is an almost essential component when
processing data using this methodology.
In this regard it is implicit,
of course, that the analyzer being used is calibrated to sample along a
sufficiently fine energy grid so as to allow the difference between signal
and noise characteristics to be discerned in the frequency domain.
The
results shown in Section 1 were obtained using high gain settings over a
limited energy band.
Since analyzers with larger storage capacity are be-
coming available, it should become a routine matter to use high gain cali-
brations over broad energy ranges, and increasingly successful separation
of signal and noise in Ge(Li) data should become possible.
In Sections 7 and 8 the extremely unstable nature of the resolution
improvement problem is clearly revealed and some promising developments are
in evidence.
The family of regularization methods has been found to offer
some promise for controlling the degrading effects of noise, roundoff
error, etc., and, in particular, the COR method outlined in Section 6
appears to be especially capable of resisting the oscillatory effects of
noise.
In all of these methods there appear in the basic formulation one
11
-------
or more parameters which effect the compromise between the considerations
of residuals and size or smoothness.
For the present time the inclusion
of these parameters appears to be indispensable but one would, of course,
prefer that such need not be chosen in an interactive fashion.
There is
some encouraging evidence in the literature [32-34] which suggests that
one may eventually be able to devise mathematical routines which allow a
good selection of the smoothing parameters from the data of the problem.
This would have the most desirable effect of eliminating human interaction
with the program and permit routine production-like application of the
methods to take place.
The imposition of nonnegativity constraints which was further required
in order to stabilize the resolution improvement calculations in Section 8
seems to be completely necessary at the present time.
Without the use of
such a priori information none of the methods considered thus far is capable
of yielding significant levels of resolution improvement of the spectral
analyzer data without incurring degradation (due to noise effects) un-
acceptable for the applications eventually desired.
Happily, this causes
only minor inconvenience and small additional expense in conducting compu-
tations.
It should be noted that the principle of imposing nonnegativity
constraints
in order to achieve a higher level of numerical stability
should find beneficial applications in a wide variety of physical measure-
ment and data transmission problems in which positive-valued intensity
data are the subject of concern.
Larger scale calculations with these
methods will be conducted during the coming year of extension of this
contract endeavor.
12
-------
RECOMMENDATIONS
An overall analysis of the results of work conducted under this con-
tract effort suggests certain general recommendations concerning the
direction of work dealing with improvement of techniques for processing of
gamma radiation data.
First of all, continued progress in the power of the
Ge(Li) analyzers is of central concern in regard to the problem of enhancing
gamma data since increased analyzer power is reflected in an increased
ability to define the photopeaks distinctly and thereby allow successful
separation of signal and noise (through methods such as those described ir,
Section I).
A somewhat paradoxical consequence of these observations,
which has not yet been adequately considered in the literature of nuclear
detection, is that although the achievement of very high levels of system
resolution typically leads one to anticipate that a photopeak will be re-
corded in only a few (perhaps one) data channels, the requirement for compu-
tational tractability necessitates that photopeaks be instead recorded in a
significantly larger number of channels.
Similar comments apply to the
resolution improvement methods in Sections 5 and 6.
Therefore, it is im-
portant to include the instrump.ntation capability, both in terms of system
resolution (keV at FWHM) and system gain (keV per channel), explicitly when
discussing any path to progress in gamma data reduction techniques.
In regard to mathematical and computational aspects of continued work,
it would seem that the constrained regularization methods represent the
most versatile and promising family of techniques for resolution improve-
ment of the gamma spectral data amonq those methods considered thus far.
Thus it is recommended that the methods of Sections 5 and 6 be tested on a
larger scale using, in particular, more complex gamma spectra as the sub-
ject data.
In order to make such increased level of complexity of test
13
-------
computations practicable it will be necessary to make a software conversion
of all procedures from the M-LAB language to a less versatile but more in-
expensive language such as FORTRAN.
This unfortunate situation is due to
the small size of the matrix calculations permitted by M-LAB and, espe-
cially, due to the exorbitantly high cost of using this program on the ADP,
Inc. network.
In the course of this work the capabilities and limitations
of the data spreading and constrained regularization techniques should be
investigated more thoroughly and alternate methods of stabilizing the calcu.
lations considered if necessary.
Throughout this contract work only single step and finite iterative
methods (singular value analysis) have been considered.
This is in accord-
ance with the major emphasis of the literature of the field of resolution
improvement problems in general.
There are, however, certain iterative
procedures which could be advantageous, and it is recommended that these
be given some consideration.
In this case it is anticipated that M-LAB
computational costs may prove to be not unduly excessive and thus this
path of quick feasibility study of iterative schemes is recommended.
For many types of early stage software development processes the M-LAB
program can be extremely beneficial.
The matrix capabilities are so con-
venient and facility with the language so easy to acquire that initial
feasibility studies concerning linear methods can often be conducted ex-
tremely quickly.
It is very unfortunate that the use of M-LAB is currently
confined to the PDP-10 family of computers and it is recommended that
support be given to the structuring of an M-LAB program in a more portable
compiler.
Finally, due to the basically visual nature of the enhancement task
14
-------
addressed in this contract work, graphics capability is a simple necessity.
As mentioned previously, the lack of high speed graphics capability at the
laboratory site has constituted a severe limitation and major obstacle of
progress of the contract work.
Therefore it is recommended that the
Eastern Environmental Radiation Facility be fully supported in its attempt
to acquire quality high speed graphics with large color screen capability
and the associated software services.
This equipment would also, of course,
prove extremely beneficial in countless other aspects of the work of the
laboratory besides those related to gamma data reduction.
15
-------
SECTION 1.
NOISE REDUCTION VIA FREQUENCY DOMAIN ANALYSIS
In this report, the term "linear filter" is used to refer to a process
L which operates upon a given input signal f to produce the output signal
g = Lf in accordance with the following model:
00
~ l~) = ~ "" le. t E. ) ~ l E) A t.
-QI
.
,
(1.1)
here the kernel function h is characteristic of the particular linear
filter in question.
In this section the additional hypothesis that h is a
difference kernel will be stipulated; thus h(e,E) = k(e-E) for some
function k.
Under this latter hypothesis (1.1) becomes an equation of con-
volution type and one often writes g = k*f to represent this situation.
If
k(x) = 0 for all x < 0 then the process L is called "realizable" since then
the output signal level g(e) (at a given energy level e) is dependent only
on the behavior of the input signal at energy levels less than e.
Although
this assumption of realizability is almost always made in usual filtering
work concerned with time-variant signals, it is not applicable to the pre-
sent endeavor and, in fact, none of the filters which are of interest here
are such realizable ones.
General background on linear filters can be
found in [2, 17, 18]; work related to that in this section has appeared in
[3, 16, 33].
The Fourier transformation concept holds the key to the filter design
process which has proved fruitful in this work.
If g = k*f and capitol
letters are used to denote the respective Fourier transforms then one has
G = KF [2].
Specifically. if
~ l~) -::
\-.. j le.\ u~ \. -1 'hoc H) c\e.
16
-------
then
"lc;,'):: \(lE) FLf)
,
(1.2)
thus the unwieldy convolution process originally indicated can be replaced
by the much more tractable multiplication operation, as indicated, by pass-
ing to the "frequency" domain via the Fourier transform.
It is assumed that the signal f, representing the output from a spec-
tral analyzer, is the sum f = d + n where d is the desired "pure" spectral
information and n is the noise superposed on d due to random variations in
the analyzer behavior, etc.
Again following the above convention with
respect to capitol letters and Fourier transforms one obtains F = 0 + N
(due to the linearity of the transform).
It is to be expected that F(E) is
approximately equal to N(E) for larger values of E; this is due to the
expectation that D(E) will be negligibly small for large E due to the
highly correlated nature of d.
This expectation forms the basis of the
following procedure for increasing the signal-to-noise ratio.
The impulse response k of the filter is constructed by means of the
following considerations relative to its transform K:
(1)
(2)
(3)
For small values of E, K(E) = 1
For large values of E, K(E) = 0
Beginning at a specified "cutoff frequency" E , the values of
c
K change from 1 to 0 along some very smooth curve.
(This last consideration is motivated by the desire to avoid introducing
oscillation related to the Gibbs phenomenon in the graph of the output
signal g.)
Then forming the product G = KF (which is equivalent to per-
forming the convolution 9 = k*f) will sharply attenuate the values F(E)
17
-------
for E significantly larger than E while retaining those values F(E) corre-
c
sponding to E < EC. By interactive choice of EC and the transition para-
meter described below, one hopes in applying this method to increase the
signal-to-noise ratio of f while suffering no significant loss of resoluti9n
(This enhanced version of f is the g above.)
Details of the Filter Construction and Application
Given the data signal f (as recorded by the spectral analyzer) one
first chooses evenly spaced samples f. thereby determining the sample
J
-+
vector f = [fo' f1' ..., fN-1]. Then, in principle, the discrete Fourier
-+
transform C'DFT") F = [Fo' Fl' ..., FN-1] defined by
r. :=.
~
£ ~\t t-x~ ll1\ \.'1 \t IN)
\t
is determined.
(This is one of a number of trivially different forms of
the OFT which could be used.)
Due to the large computational times in-
curred in the indicated formulation of the OFT (for large N) the actual
computation is carried out by means of a fast Fourier transform ("FFT")
[26]; this is true of all transforms and inverse transforms referred to in
this report and no further explicit mention of this point will be made.
The numbers F. are complex but the object of interest in the present
J
context is the plot of the quantities IF.I (the complex modulus) versus
J
the index (I' frequency number") j; in some instances it may be more appro-
priate to plot the logarithm of the modulus in lieu of the modulus due to
the large range of values typically of interest.
If the hypothesis that
f is the sum of a highly correlated signal d and a noise signal n is
satisfied then one will typically observe a transition from the higher
18
-------
power associated with the smaller frequency numbers to lower levels of
power at the larger frequency values.
This region of transition is the
likely location of the optimum values of the cutoff parameter [18]; the
specific choice finally settled upon for a given set of data is accom-
plished interactively.
The final filter design consideration involves the specific manner tn
which the transition of the values of K from 1 to 0 takes place.
The
method which was finally selected for use calls for this transition to
occur along a cumulative normal distribution curve with standard deviation
given by the parameter cr.
Thus, given the choices of the parameters EC
and cr the (piecewise continuous) description of K is given by
\ ~~+"(/
(1.3)
1. - ~ f ~ e LX ~ ~ - ~ " <.. 'j. - e < - ." f)' 1 ~.
-w
(.h l .
In the program written to accomplish the filtering process the discrete
sampling of K is accomplished using a "look-up table" to obtain aporoximate
values of the cumulative normal distribution curve ordinates.
Finally, the OFT G of the enhanced version 9 of f is obtained by
forming the products G. = K.F., and 9 is then determined by means of the
J J J
inverse OFT as
~i ~ * f (,,~ "tl-l~"i~ IN) .
+ +
Under the indicated hypothesis, g will comprise an enhancement of f having
increased signal-to-noise ratio.
The optimum values of E and cr are found
c
19
-------
interactively; one finds that the results obtained are normally quite in-
sensitive to small variations in these parameters.
This methodology was first applied to the 1024 channels of light
intensity data indicated in Figure 1.1.
(This data was kindly provided by
W. F. Ranson; it resulted from work conducted in the preparation of [25]).
The behavior of the 513 values of the discrete Fourier transform F of this
data is indicated by Figure 1.2.
In this figure one may observe the ex-
pected tendency for the modulus (which is associated with the power of the
signal) to steadily decrease with increasing frequency number until a
"leveling offll is eventually experienced in the vicinity of approximately
frequency number 15 to 25.
This behavior suggests that the cutoff frequency
€ for the filter to be applied to this data be chosen in the region from
c
15 to 25.
After some experimentation the value € = 15 was chosen as
c
appropriate.
(The results finally obtained are quite insensitive to small
changes in € .)
c
with which the transition in the values of K occurs was selected inter-
Similarly the standard deviation 0 specifying the sharpness
actively.
The value 0 = 2 was found to be adequate but, again, the results
are very insensitive to perturbations in the values of 0.
The transfer
function K corresponding to this choice of the parameters €
c
and 0 is in-
dicated graphically in Figure 1.3, and the modulus of the corresponding
result G from (1.2) is shown by Figure 1.4.
This G represents the result
of the filtering process in the transform domain, and its inverse DFT g
(Figure 1.5) thus comprises the end product of the smoothing process.
The noise reduction technique was also applied to a noisy gamma
photopeak; the raw data is shown in Figure 1.6.
In this case there were
60 channels of data and thus 31 Fourier transform values (figures detailing
20
-------
the intermediate steps have not been included in this instance).
The same
sort of consideration of the modulus of the transform as was previously
employed suggests that cutoff frequency number be taken at channel 6.
In
this case a simple truncation of the subsequent values of the transform
was sufficient to yield very acceptable results.
The results of this
smoothing is shown in Figure 1.7.
The results achieved in these two
examples (and in other applications which were not included here due to
their entirely similar nature) are especially good because of several
factors.
First of all, the modeling assumption that the raw data signal f
should consist of a very smooth deterministic component d with a much smaller
amplitude, highly erratic noise signal n superposed is especially well satis-
fied in these example problems as is rather clearly revealed by the plots
of raw data (Figures 1.1 and 1.6).
Secondly. the signal d (which may be
assumed to be accurately approximated by the result g indicated in Figure
1.5 and 1.7, respectively) is one which lends itself especially nicely to
approximation by trigonometric polynomials of low degree, and thus the
Fourier approach is exceedingly appropriate.
In any case, the noise re-
duction technique outlined in this section has been amply demonstrated to
be a powerful tool for certain types of data, and this fact has been
recognized by the appearance of the article [6].
21
-------
>-
fo-
-
C/'1
Z
~
fo-
Z
-
f0-
X
U
-
...J
A
~ I! ,
~~ ! ~ h
o J r~\ 'I Ii'
V') i I \ J \
I J \ I \ \ ~~.
~ j)K~JJ ~ v k'~1~~JJ~~
g~ lit' ~ . "
JI ~. 'll
~ f ~~
I A :
o I J~
-V
o . I
o 10 20
o~
00 !
I
I
I
01
r-
o
\0
Figure 1.1
J] I I I
30 40 50 60 70 80
CHANNEL NUMBER ~X101)
I
90
,
100
Light intensity data ([rom [ 25 ] ).
22
-------
--
>.
....
.- r-
c;t)
c
~
....
c ..0:
.... i
.c
b()
.-
- VI
"-
0
E
a-. ~
~
c;t)
c r"j(
to;!
a-.
.... I
a-. : . ~f~~~'#i~\Wi~~~~I~f~
.~
::I
0
u..
--
r./')
~ I I
-<: I :
0 .0 ~ i
b() I
0 I
....:I ' I
-~ ,
'0 10 20 30 40 50 60 70 80
frequency number (X 10')
Figure 1.2
Modulus of the transform of data from Figure 1.1
23
-------
tV
"'0
::I 0
.....
.- "":1 "\
c:
b{) \
~ \
E: 00 \
.... 0
aJ \
... \
- \0
1i: 0 \
\
~ \
01 \
NI \
. J
01
01
'- r
0 10 20 30 40 50 60 70 80 90 100
Frequency number
Figure 1.3
Filter function chosen for the enhancement
24
-------
--
;;..,
.....
r;/j
c:
(1)
.....
c:
.....
..c
QI)
.-
-
'-
o
8
I-
o
'-
r;/j
c:
Ct!
I-
....
I-
(1)
'i:
:::
o
~
--
C/"J
~
<::
o
QI)
o
..J
t-l
\0
Ir\
~
f"") "I'r,
~
N j~1
I .
- J ~i
i i
i
,
o i
I
,
-- J :
10
Figure 1.4
10
20
30 40 50 60 70
Frequency number eX 101)
80
Attenuation of transform shown in Figure 1.2.
25
-------
o
00
0
r---
>- 0
\0
r-
- 01
CI"J
Z V')-I
I
~ I
r- I
z 01
- ~1
r- I
== oj
0 r-"\I
-
..J 01
Nl
i
=:l
01
0 10
(\ r\ r\
I .\ i '\ I \
I I \ I \
I J V \ f\
~ V~
\\
\vI
Figure 1.5
20
I 1 I .
30 40 50 60 70 80
CHANNEL NUMBER (Xl0I)
9u
100
Result of enhancement of data in Figure 1.1.
26
-------
o
o
o
Lf)
-
o
o
.
Lf)
C\J
-
o
~o
Ie
(;)0
.....,-
W
I
0
o
.
Lf)
r--
o
o
.
o
Lf)
o
o
Lf)
C\J
o
o
~.oo
45.00
52.50
22.50
30.00 3 .50
CHANNEL
60.00
7.50
15.00
Figure 1.6
Raw photopeak data from a radium energy line.
27
-------
o
o
.
o
Lf)
-
o
o
Lf)
C\J
o
1-0
Io
t:Jo
.....,-
L.L.I
:r:
o
o
.
Lf)
1'-
o
o
o
Lf)
o
o
.
Lf)
C\J
o
o
.
~.oo
7.50
15.00
22.50
Figure 1.7
30.00 37.50
CHANNEL
45.00
52.50
60.00
Result of enhancemellt of data in Figure 1.6.
28
-------
SEcn ON 2.
GENERAL CONSIDERATIONS IN RESOLUTION IMPROVEMENT
The operational characteristics of a radiation spectrometer are
usually modeled by the equation*
v
~ l~) "C ~ '" l ~ , t) \ l~) d t.
~
(2. 1)
where g(e) is the recorded intensity of the radiation at energy level e,
f(E) is similarly representative of the incident radiation, and h is a
function determined by the particular spectrometer (including its settings
and calibration, etc.).
Often one makes the further assumption that h is
a difference kernel (i.e., h(e,E) = k(e-E) for some function k); it is well-
known, however, that this assumption is likely to be a good approximation
over only a limited energy range.
With this assumption the equation be-
comes one of convolution type and one may often write g = k*f (or some
similar notation) to denote this convolution process.
Equations such as (2.1) are categorized as "Fredholm integral
equations of the first kind."
There is a large amount of literature devoted
to the problem of generating reliable approximate solutions of such equa-
tions and, in this literature, the many perversities of this problem are
evident.
The essential point, however, is easily seen by examining the
following considerations [22].
For any integrable kernel h it can be shown that the values of the
function g given by
m
*In practical applications there is typically no loss of generality in
specifying a finite upper limit U for the indicated integral.
29
-------
\1
~\'W\ l~) ': ~~ '" lt~ E) siV\ l~ ~) d ~
approach 0 as m becomes large.
This means~ assuming that the map of f onto
g given by (2.1) is one-to-one~ that replacing the data g by g + gm (for
some large m) yields an equation whose solution differs from that of (2.1)
by the signal sin(mE).
Th"
~n extremely small pertubation g in the data
m
of the problem gives risl Lo a very significant oscillatory disturbance in
the (theoretically correct) solution.
In practice~ of course~ one will never have data g which agrees
exactly with that which the analyzer model would describe for the given
incident signal f.
Furthermore~ even if such "exact" data were available~
computer roundoff errors will be incurred in any practical calculation.
Even though these errors will generally be small ones~ the result of their
presence will often be a very significant (usually oscillatory) disturbance
of the true solution f.
This difficulty causes any "straightforward"
attempt to achieve an approximate solution to (2.1) to be generally unsatis-
factory. and thus methods which are (at first glance) quite unintuitive are
entirely predominant in this field.
The next several sections will discuss
the various methods which have been applied to the spectral resolution
improvement problem of (2.1).
The following notation and specific problem
statement will be relevant to a reading of these sections.
Consider the linear operator % defined by
l~ ~) l{)
":
v
~ '-'l~,E) ~lE:) At
Q
(2.2)
,
30
-------
then the formal problem associated with (2.1) can be stated as
Problem P*:
Given g and h, solve the operator equation ~f = g for f.
On the other hand, any practical computational version of this problem
would have to be phrased more nearly as
Problem P:
Given noisy samples from g and samples from h, estimate
sample values of f given that ~f = g.
In the course of this contract effort a great deal of literature re-
lated to the solution of problems similar to that given in (2.1) was sur-
veyed and the many techniques for solution of problems such as Problem P
presented in this literature were considered vis-a-vis the enhancement
problem for gamma radiation spectral data.
In the following sections the
main types of solution procedures are described and implemented on proto-
typal example problems (in Section 7 and Section 8).
In particular, the
necessity for some type of stabilizing procedure and even a priori assump-
tions (based on physical interpretation of the problem) will become evident
31
-------
SECTION 3.
QUADRATURE METHODS AND SINGULAR VALUE ANALYSIS
The most immediate (and least successful) method for Problem P involves
choosing quadrature weights w. and replacing (2.1) by
J
~\. '::
N
Z ~. ~.. \.
~~\ \ "\ ~
,
1 ~
"
~
t-\
(3.1 )
where g. is an approximate observation of g(e.), f. = f(E.) and h.. =
1 1 J J lJ
h(e.,E.). (Here the e. and E. are the sample points.) This can be written
1 J 1 J
-+ -+
g = Af
(3.2)
where A = [A..] is the MxN matrix determined by A.. = w.h.., M is the num-
lJ lJ J lJ
ber of samples of g observed and N is the number of quadrature nodes E..
J
In this context M ~ N will always hold, and it will be further assumed that
A has rank N although more will be said concerning this assumption later.
In all sections following this one the assumption that M = N will be made
for simplicity.
-+ -+
If M = N then one has f = A-1g as the formal solution.
More genera 11y.
when M : N the usual least squares criterion would be the natural sense in
-+ -+
which to consider the solution of the possibly overdetermined system Af = g.
In either case, one may use the concept of the singular value decomposition
[19] of a matrix to computational advantage.
Under this decomposition A be-
comes factored into the form UDV~ where U is an MxN matrix with orthonormal
columns, D is an NxN diagonal matrix whose diagonal entries are the (posi-
tive) singular values of A, and V is an NxN orthogonal matrix; here the
prime (~) symbol is used to denote the transposition operation.
The normal
-+ -+ -+-+
equations for the system Af = g are A~Af = A~g.
The latter is a non-
32
-------
singular system of N equations in the indeterminants fl' ..., fN but the
coefficient matrix A~A is expected to be very ill-conditioned and thus the
solution of this system must be approached with great care.
Using the
-+
singular value decomposition A = UOV~ and changing the variables f to f =
V~f and the data 9 to 9 = U~g, one obtains the diagonal system Of = g.
A -+ -+
This can be easily solved for f; f is then available from f = Vf.
Reliable solutions are generally not available from this procedure due
to the ill-conditionedness of A.
The singular values o. of A (which comprise
1
the diagonal entries of 0 - we assume, without loss of generality that they
are in decreasing order) will typically approach zero rapidly as i approaches
N.
Since 9 consists of noisy samples, one will find that the components of
A
g are sums of the form d. + n. where d. is representative of the pure spec-
1 1 1
tral information and n. results from the noise. One hopes that the d.
1 1
approach zero faster than the o. but the same is certainly not to be ex-
1
pected of the n..
1
When the respective components of f are computed from
the diagonal system Of = g the quantities d./o. + n./o. result. For larger
1 1 1 1
values of i the ratio n./o. will tend to increase wildly and hence the
1 1
noise component of the data will dominate the solution.
This analysis of the situation suggests also the following possible
remedy.
If V. is the j th column of V and f has components c. as described,
J J
least squares solution to Af = 9 can be expressed as f =
thrn the
N -+
I c.v..
j=l J J
As mentioned above, it is anticipated that the latter summands
in this expression will be dominated by the effects of the noise
-+
component of the data g and thus one may wish to consider deleting some of
these.
The term "singular value analysis" refers to the
process of
k -+
= I C.v.
j=l J J,
inter-
-+
actively considering the behavior of the partial sums fk
33
-------
k = 1, 2, ... N.
In this regard the positive integer k can be considered
as a "smoothing parameter" with smaller values of k giving high signal-to-
noise ratio but low resolution level while the larger values of k yield the
reverse compromise situation.
It should be observed that this process
amounts to the truncation of the last N-k singular values of A thereby
producing a computational version of A having rank k.
This technique was applied to several relevant resolution enhancement
problems during the course of this contract effort.
The program of Golub
and Reinsch [12] was used to perform the singular value decompositions of
the various matrices as required while the software otherwise necessary
was implemented on the M-LAB program of the ADP, Inc. Network.
Typical re-
sults obtained by this means are indicated in detail in Sections 7 and 8.
34
-------
SECTION 4.
GENERAL PHILOSOPHY OF REGULARIZATION METHODS
There are a number of methods in the literature which may be loosely
identified as "regularization methods."
Roughly speaking, each of these
methods attempts to provide a compromise between the residual in (2.1)
on the one hand, and the size of some norm of the solution f, on the other.
(This norm may be a reflection of the size of f or one of f's derivatives,
or some combination of these.)
The philosophy of regularization methods is fairly well revealed by
the statement of the following problem.
Pro b 1 em R:
Given the operator~, the data g (possibly contaminated
by noise), the norm 11.11 R' and the number A :: 0,
determine f = fA so as to minimize the expression
2 2
Ilg- fll + A211fll . Here the expression 11.11
2 R 2
to the L (RMS) norm given by
2
refers
\\ ~ \\). ':
'h.
l ~~ ~\~) A E 1.
In this class of methods the "smoothing parameter" A is often chosen
interactively (as was true for the singular value analysis parameter k of
the previous section).
However, there has also been some fairly recent
work [32] in which the choice of this regularization parameter is made by
means of specific considerations concerning the data 9 and assumptions
relative to the characteristics of the noise which is contained within
this data signal.
Due to the fact that the noise models considered thus
far in this type of work are not appropriate for the noise encountered in
35
-------
the spectral data from the analyzers of interest in this contract endeavor
and to the preliminary nature of the work reported herein, these methods
for the choice of the "optimum" value of A were not actively considered in
the computations which were carried out.
The interested reader is referre~
to [32, 33, 34] and the references contained therein.
In Section 8 it will be demonstrated that the application of the a
priori knowledge of the nonnegativity of the functions being considered can
have a most powerful effect on the quality of the solutions obtained.
Thus
one may find it beneficial to modify the statement of Problem R to require
the minimization of the indicated form among functions f which have no
negative values.
In such a case as this the problem of attempting to
select the "optimum" value of the regularization parameter is one of active
current research interest.
36
-------
Section 5.
Discrete Implementations of Regularization
In this report the term IIdiscrete regularizationll (the terminology of
this field is not at all standardized) will be used to refer to a form of
regularization in which the quantities in the statement of Problem R
(Section 4) are all calculated using discrete operations.
In this section
several methods based on the discrete regularization concept will be formu-
lated and their application to sample problems studied.
Euclidean Regularization
The term IIEuclidean Regularizationll (IIERII) will be applied to methods
based on the following idea.
-+ -+
g = Af of (2.1) (which one obtains by applying a quadrature formula to the
Considering again the discretized version
indicated integral) one may consider the problem:
given A~O, minimize the
expression
-+-+2 -+
Ilg-AfiIE + A211fll~
( 5. 1)
-+ -+
by properly selecting f = fA. In this (discrete) context the norm II-liE
is that of the Euclidean length in N-space, i.e., IlfilE is defined as
\?: ~ 'i r"
~
This problem may be solved rather efficiently [31] by realizing that
-+ -+ 2 -t: 2
the terms Ilg-AfiIE and A211TIIE can be thought of as the total square
-+ -+
residuals in the linear systems Af = 9 and AIf = 0, respectively. It
follows, therefore, that the minimization problem stated above is identical
with that of determining the least squares solution to the linear system
L\11 = \\\
37
-------
and hence the solution, in principle, is available from the associated normal
equations
~A,d\ \~l\ ~
=
\~,~ll \,'
which can be written as
-+ -+
[A~A + A2I]f = A~g .
(5.2)
It is advantageous to write A (again as in Section 3) in the form
Then A~A + A2I = V[D2 + A2I]V~, and 50 if one sets f = Vf and
UDV~ .
-+ ~
g = Ug then (5.2) can be written
[02 + A2I]f = Dg .
(5.3)
Due to the diagonal character of this latter system the components of fare
immediately available by performing the obvious divisions; then one easily
-+ -+
calculates f = fA by the multiplication of V and f.
calculation of fA (for the various values of A which
In this way repeated
one may want to consi-
der in a typical interactive situation) can be very quickly and inexpensive-
ly generated.
Differential Euclidean Regularization
Other forms in which discrete implementations of regularization may be
advantageously cast involve formulations dependent upon the difference
operator ~ given by (~f). = f"+l - f..
1 1 1
It should be observed that ~ is a
linear transformation from N space onto (N-l) space and, thus, that the
-+
second order operator ~2 given by ~2f = ~(~f) maps N space linearly onto
N-2 space.
~ and ~2 are, of course, often used in constructing discrete
approximations to the first and second derivatives, respectively, of a
function whose samples f. at equally spaced grid points have been determined.
1
38
-------
Due to the widely observed fact that approximate solutions to (2.1)
are often subject to pathological oscillation, it may be helpful to regu-
larize the solution of the discrete approximation AT = 9 in one of the
following ways.
First, one may consider the optimization problem:
given
A, g and 0, A~O, minimize the expression
119 - Afll~ + A211fll~ + (AO)2116fll~
(5.4)
One could view this as the method which results from the replacement of the
-+
second Euclidean norm in (5.1) by the norm 11.110 given by Ilfllo =
[llfll~ + 02116fll~]\ Thus, this "First Order Differential Euclidean
Regularization Method" ("DERl") is concerned with the minimization of
-+ -+ 2 -+ 2 -+
Ilg - AfllE + A211fllo (given g, A and A, o~O) and hence results from a
small, formal (although important) alteration of (5.1). The same sort of
consideration as was given to the ER optimization problem now shows that
-+ -+
the DERI solution to Af = g can be realized as the least squares solution
to the system
\A'~ \ ~ = l \ ]
(5.5)
where, in this context, "6" now denotes the matrix of the linear transforma-
tion 6 relative to the standard bases of the Euclidean spaces of dimension
Nand N-1, respectively.
(Specifically. then the jth column of the matrix
-+
6 is the image 6e. of the jth column of the NxN identity matrix I under
J
the transformation 6.)
To solve tris problem one may write the singular value decomposition
of the 3N-1 by N coefficient matrix of the system (5.5) as UDV~.
In this
case, the three factors are functions of the parameters A and 0; U is 3N-1
by N with orthonormal columns while D and V are both N by N.
-+
If f is
formally set to be Vf and g is defined as the product of U~ and the vector
39
-------
on the right hand side of (5.5) then the normal equations for (5.5) can be
transformed to the diagonal system Of = g.
(Here U~U yields the N by N
identity matrix due to the fact that the columns of U are orthonormal.)
-+ -+
Then the coordinates of f can be easily found by division and f = fA,o
A
generated as the product Vf.
Another type of differential Euclidean regularization which can some-
times be used to advantage in suppressing oscillation in the solution of
(2.1) employs consideration of the growth of the second order difference
-+
~2f in the formulation of the solution. This "Second Order Differential
Euclidean Regularization fv1ethod" calls for replacing
A1 = 9 with that of minimizing the expression
-+ -+2 -+2 -+2
Ilg - AfllE + A211fliE + (AE)2 11~2f11E
the problem of solving
(5.6)
In a manner completely analogous to that used following (5.4), one could
consider this IDER2" method to consist of the minimization of the quantity
-+ -+2 -+ 2 -+
Ilg - AfllE + [AllfIIE] where the norm II-liE is defined by IlfilE =
2 2 k:
[llfllE + E211~2fIIE]2.
-+ -+E -+ -+ -+
The solution f = fA to Af = g (given A, g and A, E~O) is, of course,
the least squares solution to the system
\:1 \ ~ "= \~ \
}.£K ()
(5.7)
where, here, 62 denotes the N-2 by N matrix whose columns are the respective
images of the second order difference operator on the standard basis vectors
for N space.
This is a system of 3N-2 equations in the indeterminant coor-
40
-------
+
dinates of f and, using exactly the same sort of procedures as was followed
in deriving the OERI solution, one can transform the normal equations for
(5.7) to a diagonal form such as Of = g.
+ +s
Then the solution f = fA is avail-
able as the product Vf.
(Here, of course, the singular value decomposition
of interest is that of the coefficient matrix of (5.7) into the form UOV~
where, in this instance, U is a 3N-2 by N matrix with orthonormal columns.)
Obviously, many other variations along the lines which have been illus-
trated here could be similarly formulated and implemented, but there does
not seem to be any justification in the context of the present study for
persisting further in this.
It seems intuitively clear that the differ-
ential Euclidean regularization methods yield a significant generalization
of Euclidean regularization and that in some types of problems this added
generality may bring welcome improvement of quality to the solution.
It
should also be clear that this increased capability has not been achieved
without the cost of additional computational inconvenience.
The singular
value decompositions must now be performed on matrices which have almost
three times the number of entries as previously; also, there are now two
regularization parameters to consider instead of the previous one.
Further-
more, the differential methods (as implemented) require that a singular
value decomposition be performed for each choice of the smoothing parameters
whereas the Euclidean method needs only one such decomposition (for a given
coefficient matrix A).
Therefore the cost incurred in typical interactive
computations with these three discrete regularization methods is normally
significantly lower when using the first of them (Euclidean regularization),
although it would be possible to reduce the cost of the differential methods
by employing constructions due to Elden or van Loan [31] (as mentioned in
41
-------
the Summary).
Finally. it should be observed that any of these regularizations can
be implemented with one or more constraints imposed.
In Section 8 it is
demonstrated that Euclidean regularization with nonnegativity constraints
can be a very attractive alternative.
42
-------
Section 6.
Continuous/Discrete Regularization
The term "Continuous/Discrete Regularization Method" ("CDR") will be
used to refer to a practical computational method [32] based on a theoreti-
cal formulation which represents a compromise between the exact solution of
(2.1) (which would be normally totally impracticable) and the discrete
regularization methods.
The basic idea of this method and the approach to
its solution are indicated in the following.
If ~ denotes the linear integral operator defined in (2.2) then given
(noisy) approximate observations g. ~ g(e.). 1 < i < N and A > 0 one con-
11-
siders the problem of minimizing the expression
...L
~
;: \\'\\~ )le~) - ~. \'l
~..\ ~
-\-
\1
). ~ ~\,,\ H.
o
(6.1)
-+
by appropriately selecting the function f = fA, If r denotes the vector of
residuals given by r. = (~f)(e.) - g. and 11.11 denotes the usual L2
J J J 2
functional norm (defined in Section 4), then the problem becomes that of
-+ 2 2
minimizing II rilE + NAllfl12 and hence the mixture of discrete and continuous
considerations becomes apparent. The solution f of this problem can be
determined (in closed form) by considerations which will now be briefly
outlined.
(The kernel h of (2.2) is assumed to satisfy the hypothesis
stipulated in [32].)
If the functions hk are defined by hk(E) = h(ek,E)
then the solution f = fA can be written as the linear combination
n
f = L ckhk'
k=l
scalars c1' ..., cN which minimize the function ~ given by the expression
Therefore the problem becomes that of determining the
43
-------
N
~(c1' ..., cN) = I [
j=l
N 2 2
(I Ckhk)(e.)-g.] + NAllIckhkl1
k= 1 J J k 2.
The quantities ck of interest can, in principle, be determined from the
extremal conditions requiring that all the partial derivatives of ~ be
zero when evaluated at the points c1' ..., cN.
In determining the solution to these equations the following notation
is employed.
-+ -+
g and c denote the N vectors whose components are the
observations g. and the coefficients c. being sought, respectively. Also,
1 1
the N by N matrix M is that given by
v
M.. = (' h (e. , E) h (e. , E) dE.
1 J Jo 1 J
-+ -+ -+ -+
Then the extremal conditions are reflected by the system M(Mc - g + NAC) = 0
(6.2)
and, using the assumption that M is nonsingular, this becomes
-+ -+
[M + NAI]C = g.
(6.3)
Solutions of this system of equations for the ck then yields the COR
solution f = fA to (2.1) in the form f = I ckhk. This is in marked con-
trast to the results generated in Section 5, of course, since one is now
determining not approximate samples of the solution f, but, instead, a
closed form expression approximating f.
A significant computational improvement over the obvious least squares
-+ +-+
solution of (6.3) as c = [M + N A I] g (where the plus sign denotes the
pseudoinverse) can be achieved in the following manner.
Since the matrix
M is symmetric, the spectral theorem of linear algebra guarantees the
existence of a decomposition of the form VOV~ where V lS an orthogonal
matrix and 0 is a diagonal matrix (the diagonal entries of 0 are the
eigenvalues of M). Using this decomposition, (6.3) can be written as
-+ -+ -+ A -+ A
V[D+NAI]V~c = g and so if one lets c = Vc and g = Vg then the diagonal
system
44
-------
[D+N\I]c = g
(6.4)
results.
A
The entries of c are then immediately available by division, one
~
gets c = VC and then generates the solution f = fA in the form f = L ckhk-
This allows extremely efficient interactive regularization iterations since,
for a given kernel function h, the decomposition calculation (of M as VDV~)
is performed only once; then the procedure indicated by (6.3) is repeated
for the various values of A which are to be considered.
In many cases a closed form expression for the kernel function h may
not be available.
In such a case it will, of course, not be possible to
determine the entries of the matrix M in the fashion indicated above;
neither will the functions hk be available in closed form and so the repre-
sentation
f = LCkhk (6.5)
must be viewed in a different light. It is assumed that one does at least
have available the samples hit = h(eift) for 1 ~ i ~ Nand 1 ~ t ~ Nand
that a quadrature rule with weights Wt and nodes Et has been selected (as a
replacement for integrations in the variable E). Then the description of
the entries in (6.2) will be replaced by the approximate calculation
M.. ~ I Wth'th't
1J t 1 J .
Also, the availability of only the samples hit will dictate that (6.5) now
describe only a sampling of the solution f = fA, Letting ft denote the
approximation to f(Et) thus generated, (6.5) will be replaced by the
equation ft = I ckhk(Et) = I ckhkt- This sampling of f can be notationally
k ~ ~ k
improved by letting f = fA denote th~ vector of samples ft and defining a
~ ~
matrix H whose entries are Htk = hkt = h(ek' Et)' Then if c = cA is the
solution to (6.3)(in which M has now been replaced by the approximation
(6.6)
45
-------
generated as indicated in (6.6)) one has the compact description
-+ -+
fA = HCA
(6.7)
The notation chosen here emphasizes the fact that the matrices M and Hare
entirely determined by the kernel function h (and the choice of grid points
ei and Et) and, thus, the selection of various values of A to be considered
in interactive iteration merely requires that the (inexpensive) repeated
solution of (6.4) be carried out.
It follows, therefore, that the cost of
such application of COR consists preponderantly of the cost of calculating
the entries of the matrix M (H is obviously a trivial matter in comparison).
Since the determination of M by (6.2) will seldom be possible in
applied work, the cost associated with (6.6) must be considered very seri-
ously.
During this contract endeavor the calculation of such Mis of the
order of approximately 40 by 40 has required an expenditure of more than one
thousand dollars on the ADP, Inc. Network.
(Of course, this cost could be
cut tremendously by an eventual development of Fortran routines for execu-
tion on the Comnet system.)
It was possible, however, to partially solve
this economic problem in a manner which will be described next.
In this discussion it will be assumed that h is a difference kernel
given by h(e,E) = k(e-E) where e and E take on values between 0 and U;
then k(s) is defined for s between -U and U and this definition may be
extended to all other real s by setting k(s) = O.
In this situation (6.2)
becomes
For many pairs i,j the two
v
M.. = r k(e.-E)k(e.-E)dE.
lJ ~Q 1 J
functions multiplied here, when
considered as
functions of the variable E are such that the intersection of their supports
is included within the interval of integration.
(This condition will be
46
-------
referred to as the
M.. =
lJ
is appropriate and, using the
"support hypothesi S".) In such a
Q)
S k(e.-E)k(e.-E)dE
1 J
-'"
change of variable s =
case the formulation
e.-E, this can be
1
written as
M.. =
lJ
SlIIIk(S)k(S+6. .e)ds
-IJI lJ
(6.8)
where 6.. e = e. -e..
1 J J 1
Any quadrature rule could be applied to the approximate computation
of M.. in (6.8) but it is of interest here to use equally-spaced nodal
lJ
points St and equal
quadrature weights Wt = 6s. Then (6.8) is replaced by
t-1.. = 6 s I k ( s t ) k ( s t + 6 . . e) . ( 6 . 9)
1 J t 1 ,J
The important fact here is that sums of the form of (6.9) can be computed
extremely efficiently through the use of the FFT algorithms previously
mentioned in Section 1. Specifically, letting kt denote k(st) it is known
[2] that the sums ft ~ ktkt+r are the values of the inverse Fourier trans-
form of the (discrete) power spectrum of the sequence kt' t = 1, ..., N
provided that the nontrivial samples from k are augmented with a sufficient
number of zero samples.
(Additional details concerning convolution and
correlation are given in [4].)
Furthermore, the power spectrum sequence
required can also be efficiently found using the FFT since its values are
the squares of the moduli of the (complex) discrete Fourier transform values
Kt- It should also be mentioned that this indirect method of determining
the M.. satisfying the support hypothesis is not necessarily more efficient
lJ
than the direct calculation of (6.9) for smaller numbers of nodes St' but
for larger and larger problems the use of this technique provides increas-
ingly more important savings. Those entries M.. not satisfying the support
lJ
hypothesis have been computed using (6.9) although for very large numbers
47
-------
of quadrature nodes a calculation based on a convolution formulation could
provide significant savings through
the use of the FFT as above.
In principle the COR method is applicable when the solution f lies in
either the space L2 of functions which are square integrable on the interval
[O,U] or in a reproducing kernel Hilbert space with reproducing kernel
satisfying certain regularity conditions as detailed in [32].
In the next
section the results of the application of this method to a well-known ill-
posed problem of D.L. Phillips are presented.
The high quality of these
results (compared to those obtained with totally discrete methods) is
apparent.
Further, in Section 8, the application of this method to sample
resolution improvement problems for gamma spectra is considered.
In
Section 8 the problems of estimating a certain moment of the incident
spectrum f as well as the spectrum f itself are considered.
The spectrum
f is modeled as a linear combination of delta functions and thus the latter
problem does not fall within the conceptual framework of the COR method as
formulated, per se.
However, sample calculations of this type were per-
formed to gain some insight concerning the robustness of the discretizations
of the method as applied in Section 7, and good (comparative) results were
still obtained.
48
-------
Section 7.
The Problem of D.L. Phillips
A classical problem in the literature of first kind Fredholm integral
equations is that stated by Phillips in 1962 [22].
In fact, the methodo-
logy devised by Phillips for the approximate solution of these types of
equations can be considered one Jf the original regularizations methods.
The problem of interest is the following.
First one sets k and g to be
the functions defined by
L t 1 -\ (01 l'li ""'1 )
Rl",\'.I:
\J t.\Sl
,~
-'\ ~ 't ~ \ ~
(7.1 )
and
5l"" ':
l\p~)\.1~ t
~ S\Y\ l t. ). I-d
(oS C'iI "'/1 ) 1 +
C7.2)
i\
o ~ ).. ~ ~ .
~ t-..,
t.\tt. .
Then the problem of Phillips can be stated in the following way:
the function f, defined on the interval [-6,6] such that
determine
"
\ '" l-f.,s' \ts) h ':
-~
~ ( 'f.'
\ -~ ~ ~ ( L
(7.3)
where g is as given by (7.2) and where h(x,s) is defined to be k(s-x).
It
can be shown that the exact solution of this problem is the function f given
by f(s) = k(s); Figure 7.1 shows the comparison between g and f.
It should
be noted that h is a difference kernel and thus (7.3) is an equation of
convolution type; therefore a number of computational "shortcuts" are avail-
able for the solution of this problem.
49
-------
In Section 2 general qualitative remarks concerning the unstable
nature of the first kind integral equations were made, but some useful in-
sight may be afforded to a reader having no direct experience with such
problems by considering the results of a straight-forward approach to the
solution of (7.3) using the quadrature methods.
To furnish such an indica-
tion of the nature of this problem the following exercise was conducted.
The integral of (7.3) was approximated through the use of both
Simpson's rule and the rectangle rule (thereby furnishing the appropriate
discretization of this equation); in both cases there were 41 quadrature
nodes taken (equally spaced) on the interval [-6,6].
Furthermore, since
in any real problem the data are likely to be obtained from experiment, it
is not reasonable to include the exactly computed corresponding samples from
g and so these sample values (which range in size from 0 to 8) were all
rounded to one fixed-point decimal place accuracy.
The resulting discrete
problems then have the form of (3.1) and (3.2) wherein the quadrature
weights w. correspond to Simpson's rule and the rectangle rule, respectively.
J
-+
Here N = 41, and one wants to solve for the vector f.
The singular value decomposition was used (as described in Section 3)
-+
to compute the solutions f to the two respective equations (3.2); the re-
sults obtained for the sample values of f have been included as Figures 7.2
and 7.3.
These figures show that the computed "solutions" are totally
dominated by the effects of noise.
(The noise results from the finite
arithmetic of the computer, the truncation errors in the quadrature rules
and the initial data errors due to the data rounding.)
In pa rti cul a r,
while the true solution f forms a smooth cosine bell varying between 0 and
2, the computed samples are totally wildly dispersed ranging from about -88
50
-------
to 46 with Simpson's rule and -16 to 18 when using the rectangle rule.
Furthermore, it is important to note that the discrete "solutions" are
actually quite good in the usual sense of residuals.
-+
Thus the vector f
which is plotted in Figure 7.2 satisfies equation (3.1) to what would
ordinarily be considered a very good "level of accuracy" since, for each
value of i, the difference between the quantities on the two sides of (3.1)
is quite small.
The fact that this smallness of residuals does not reflect
high accuracy of the approximate solution is characteristic of ill-
conditioned linear problems.
Similar comments can be also made with regard
to the approximate solution indicated in Figure 7.3.
The singular value analysis described in Section 3 was carried out
for the problem of Phillips using both exact data and data rounded to one
fixed point decimal place; in each of these instances both Simpson's rule
and the rectangle rule were used to supply the quadrature weights w. in (3.1).
J
To lower the cost of computation the numbers of grid points were reduced to
N = 13.
The results are indicated by Figures 7.4 through 7.14; in these
figures the theoretically correct solution to the integral equation is indi-
cated by a trace of stars(*) while the thirteen computed sample values of f
are superposed by plotting zero's (0).
Whenever two points are found to
overlie one another (to the resolution level of the plotter used) the two
symbols are superposed - this convention remains in effect throughout this
report.
It is very interesting to note that Simpson's rule (which would
ordinarily be considered generally superior to the rectangle rule) gives
oscillatory results (Figures 7.4 through 7.7) which eventually (Figure 7.7)
considerably lIovershoot" the correct peak maximum.
This pathology is evident
51
-------
even with exact samples from g and the rounding of the data merely causes a
further slight deterioration of the same qualitative features (Figure 7.8
and 7.9).
On the other hand, using the rectangle rule with the same exact
and rounded data gave much more regular and accurate results (Figures 7.10
through 7.13) with the optimum approximations obtained using k = 9 and k = 7,
respectively.
(Note that the smaller value of the optimum k when higher
levels of noise are present is in accordance with the theory outlined in
Section 3.)
Previously, (when using Simpsons' rule and exact data) the
optimum value of k was also 9 but the results (Figure 7.6) were not nearly
as acceptable as with the rectangle rule (Figure 7.11).
The explanation
is that when using the latter quadrature method the resulting eigenvectors
~
v. as described in Section 3 represent discrete versions of eigenfunctions
J
shaped much more compatibly with the solution f than is true when using
Simpson's rule.
For this reason most of the emphasis in the remaining
computations to be discussed herein will be placed upon solutions in which
the rectangle rule is used for performing quadrature calculations (whenever
such are needed).
The various discrete forms of regularization discussed in Section 5
were also applied to Phillips' problem.
Both types of quadrature methods
were applied for comparison; also, both exact and rounded data were con-
sidered.
The results obtained using Euclidean regularization (Figures 7.14
through 7.20) show that the same sort of oscillatory pathology is incurred
when using Simpson's rule (Figures 7.14 through 7.17) as was previously
experienced when using singular value analysis with this quadrature method
(Figures 7.4-7.9).
Similarly. the results obtained using the rectangle
rule (Figures 7.18-7.20) are much more stable with very high quality being
52
-------
obtained with zero regularization when exact data is used (Figure 7.18) with
only a slight degradation occurring as a result of rounding of the data
(Figure 7.20).
The alternate forms of regularization were also applied to this pro-
blem and some typical results are included (Figures 7.21-7.37).
The compari-
sons here with previous results are quite interesting.
Figures 7.22 and
7.25 show that the first order differential regularization is quite capable
of suppressing the pathological Jscillation in the computed solution which
was previously invariably present (Figures 7.4-7.9 and 7.14-7.17) when using
Simpson's rule to perform the quadratures.
Unfortunately. this increased
power comes at the cost of higher computational expense (as compared to the
Euclidean regularization method, for example).
As previously mentioned, it
is possible to develop more efficient computational methods for performing
DER1 regularization, but it was not felt appropriate to expend any develop-
mental effort on this since no advantage over the ER regularization method
was observed (on these types of problems) when using the rectangle rule.
Similar comments apply to the use of the second order differential
regularization method (Figures 7.26-7.29).
The previous pathology associated
with the use of Simpson's rule is largely eliminated (Figures 7.27 and 7.29),
but the computational expense is higher and no real advantage over less
expensive techniques was noted when using the rectangle rule.
A considerable amount of effort was expended on the continuous/discrete
regularization technique; typical results are displayed in Figures 7.30
through 7.37.
It was found that the oscillation associated with the use of
Simpson's rule was very effectively resisted (Figure 7.31 and 7.33) and
that very nice results were also obtained if the rectangle rule was used
53
-------
instead (Figure 7.35 and 7.37).
Moreover, as described in Section 6, very
dramatic reductions in computational costs were achieved by various means.
Although only a slight increase in quality of results over the Euclidean
regularization method was observed when using the rectangle rule, it is
felt that the COR method deserves further considerations.
54
-------
8.91
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 0 +
o 0
8.19 + ;- Exact Data +
7.47 + "
0 +
6.75 + +
6.03 + +
o 0
5.31 + +
4.59 + +
o 0
3.87 + +
3.15 + +
o 0
2. 43 + +
+ [Exact Solution
*
1. 71 + +
o * * 0
.990 + * * +
o 0
* *
.270 + 0 0 +
.~ ~ $ * * * * * * * * * ~ $ $ $.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.1
Comparison of data and solution.
55
-------
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .T.... .+.... .+.... .+..
44.3 + +
o 0
30.9 + +
17.5 + ~ +
I AI
4.13 + \ I \
\ I I \
I I I I
I
I I
-9.27 + I ~ i
I \
I
-22.7 + I I
I
i \
-36.1 + ! I i +
\ I I
\ I \ I
-49.5 + \ +
\ I \ I \'
I \
i j I
~ I I
-62.9 + \: I I +
!
j
-76.3 + +
.+.... .+.... .+.... .+.... .+.... .T.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1. 20 3.00 4.80
Figure 7.2
Solution by quadrature method using Simpson's
rule.
56
-------
18.0
14.4
10.8
7.15
3.55
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 0 +
+
+
+
+/
. \
. \
-.500H + "
. \
-3.65
-7.25
-10.9
-14.4
+
+
+
+
!
o
i :
\ I
, I
I I
,
,
\ i
~
i
, i
'/
, I
!
II
1
\ 0
~,
I I ~\
I '
,
,
, I
I i
! '
: i
I
I
1
i
i I
I I
Ii
I'i
I I
i ,
': '
I !
: I
i 1,\
I '
, I
I '
, ,
i ~ r
'Ii j ~
, ~
i
,
if
1
+
+
+
(0\ ~
\ .
! \ +
I 6.
i
+
~
I
\
R \ I
I b I \
\\! \ II
! ~ \
\
\
\
! ,
, ,
I i
i
i
I ! I
i I ; I I' : :
I I! i
, ,
, I
1 !
, i
, I
Iii i
" ,I
+
+
+
. +. . . . .+.. . . .+. . . .. +. . . . . +. . . . .+. . . . .+.. . .. +. ....+. . . . .+.. . . .+. . . . .+.. .. . +. . . . . +. .
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
I '
I !
'I!
,I
: i
!
Ii 1:1 II
~ j 1
Figure 7.3
Solution by quadrature method using rectangle
rule.
57
II
"
\1
'I
1
-------
1. 98
1. 78
1. 58
1. 38
1.18
.975
.775
.575
.375
.175
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
'"
*
+
*
+
*
+
+
*
*
+
+
+
*
+
*
+
+
*
+
+
+
+
+
+
*
'"
+
+
*
*
.,o, * * * *
* * * *
*
* '" '" * '"
'" * * * "'.
. +. . . . .+. . . . .+. . . . . +. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. .
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.4
Singular value analysis with exact data and
Simpson's rule; k = 1
58
-------
1.97
1.74
1.50
1. 27
1.03
.800
.566
.331
.973@-1 +
-.137
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
+
+
+
+
+
+
+
.* * * * *
+
*
*
*
*
+
+
+
+
+
+
+
*
*
+
*
*
*
* * * * *.
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.5
Singular value analysis with exact data and
Simpson's rule; k = 3.
59
-------
1. 9 7
1. 76
1.55
1.34
1. 13
.918
.707
.496
.285
.735@-1 +
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
+
+
+
+
+
+
+
+
*
*
+
+
+
+
+
*
+
*
+
*
*
*
*
*
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.6
Singular value analysis with exact data and
Simpson's rule; k = 9.
60
-------
2.91
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
+
2.61 +
2.3Z +
2.02 +
1.73 +
1.43 +
1.14 +
.846 +
.552 +
.257 +
+
+
+
+
+
+
+
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.7
Singular value analysis with exact data and
Simpson's rule; k = 13.
61
-------
1.98
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
*
*
1.78 +
1. 58 +
1.38 +
1.18 +
.975 +
.775 +
.575 +
.375 +
*
.175 +
*
*
+
+
* *
+
* * +
+
+
+
+
*
+
.0 * *
.* * * * * * * * * * * * * * * * * * * *.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.8
Singular value analysis with rounded data and
Simpson's rule; k = 1.
62
-------
3.56
.+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . . +. . . . .+. . . . .+. . . . . +. . . . .+. . . . . +. . . . .+. . . . .+. .
+ +
+
..
..
+
.. ..
+
* *
+
+
..
+
..
+
3.20 +
2.84 +
2. 48 +
2.12 +
1.76 +
1.39 +
1. 03 +
.675 +
..
.315
..
+
+
..
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.9
Singular value analysis with rounded data and
Simpson's rule; k = 13.
63
-------
1.98
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
*
1. 78 +
1. 58 +
1. 38 +
1.18 +
.975 +
.775 +
.575 +
.375 +
.175 +
*
*
*
*
*
*
*
*
*
+
+
*
+
*
+
+
*
+
~ ~
~~
*
+
*
.* * * * * * * * * * * * * * * * * * * *.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.10
Singular value analysis with exact data and rec-
tangle rule; k = 1.
64
-------
2.00
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
1.80 +
1.60 +
1.39 +
1 . 19 +
.986 +
*
.782 +
.578 +
$
0375 +
*
-171 +
*
* *
+
* *
+
$ $
+
* * +
+
*
+
+
$
+
*
+
*
.$ * * $ * * $ * * $ $ * * $ * * $ * * $.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .-t.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.11
Singular value analysis with exact data and rec-
tangle rule; k = 9.
65
-------
1.99
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
*
*
1.79 +
1.58 +
1.37 +
1.17 +
.960 +
.754 +
.547 +
$
.340 +
*
.133 +
*
*
+
+
$
$
+
* * +
+
* *
+
+
$
+
*
+
*
o 0 0 0
.$ * * * * * * * * $ . * * * * * * * * $.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.12
Singular value analysis with rounded data and rec-
tangle rule; k = 7.
66
-------
2. 3 7
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
2.13 +
1. 89 +
1. 65 +
1. 41 +
1.17 +
.930 +
.690 +
.450 +
.210 +
+
+
*
+
+
+
+
+
+
*
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.13
Singular value analysis with rounded data and rec-
tangle rule; k = 13.
67
-------
2. 91
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
+
2. 61 +
2.32 +
2.02 +
1. 73 +
1. 43 +
1. 14 +
.846 +
.552 +
.257 +
+
+
+
+
+
+
+
+
J. e J. J. e.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.14
Euclidean regularization with exact data and
Simpson's rule; A = O.
68
-------
1. 97
. +. . . . .+. . . . .+. . . . .+. . . . . +. . . . .+. . . . . +. . . . .+. . . . .+. . . . .+. . . . . +. . . . . +. . . . .+. . . . .+. .
+ * +
* *
1.77 + *
1.57 +
1. 36 +
1.16 +
.955 +
.751 +
.548 +
\ *
.344 +
*
.140 +
+
+
+
+
+
+
+
+
+
. ----0--- * 0---0 ~.
.~* * * * * * * * * * * * * * * *~.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.15
Euclidean regularization with exact data and
Simpson's rule; A = .030.
69
-------
3.56
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
3.20 +
2.84 +
2.48 +
2.12 +
1.76 +
1.39 +
1.03 +
.675 +
*
.315 +
*
+
+
+
+
*
* *
+
* *
+
* *
+
+
*
+
*
*
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.16
Euclidean regularization with rounded data and
Simpson's rule; A = O.
70
-------
1.97
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
* *
1.77 + * *
1.56 +
*
1.36 +
1.15 +
.950 +
.745 +
.540 +
.335 +
*
.130 +
+
* * * * * * * * * * * ~:
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
+
+
+
+
+
+
+
+
Figure 7.17
Euclidean regularization with rounded data and
Simpson's rule; A = .030.
71
-------
1. 97
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
*
*
1. 77 +
1. 5 7 +
1.37 +
1.17 +
.975 +
*
.775 +
.575 +
*
o
.375 +
*
.175 +
*
*
*
+
+
Ii Ii
+
* * +
+
*
+
+
*
0
+
*
+
*
.$ * * . * *. * * * . * * . * * * * * *.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... ."t".... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.18
Euclidean regularization with exact data and rec-
tangle rule; A = O.
72
-------
2.37
. +. . . . . +. . . . .+. . . . . +. . . . .+. . . . . +. . . . .+. . . . .+. . . . .+. . . . .+. . . . . +. . . . . +. . . . . +. . . . . +. .
+ +
2.13 +
1.89 +
1.65 +
1.41 +
1. 17 +
.930 +
.690 +
.450 +
.210 +
+
+
+
+
+
+
+
+
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.19
Euclidean regularization with rounded data and
rectangle rule; A = O.
73
-------
1.97
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ e +
*
*
1. 77 +
1. 57 +
1. 3 7 +
1.16 +
.962 +
.760 +
.557 + 0
*
.355 +
*
.152 +
*
*
+
+
* *
0 0
+
* * +
+
* *
+
o +
*
+
*
+
* 0
o *
.. * * * * *. * * . . * * e * * * * * e.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.20
Euclidean regularization with rounded data and
rectangle rule; A = .030.
74
-------
2.20
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
1. 97 +
1.75 +
1.52 +
1.30 +
1.07 +
.847 +
.621 +
.396 +
.170 +
+
+
+
+
+
+
+
+
+
* * * * *
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.21
DERI regularization with exact data and
Simpson's rule; A = .050, 0 = 1.
75
-------
1. 97
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ e +
*
*
1.77 +
1.56 +
1.35 +
1 . 15 +
.939 +
.732 +
.525 + *
o
.318 +
*
.111 +
0 *
. * * * & * * * * * .
.0
*
*
+
+
13 e
+
* *
+
+
* *
+
* +
o
+
*
+
*
o
e * * * *
* 13 * * e.
o.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.22
DERI regularization with exact data and
Simpson's rule; A = .100, 8 = 1.
76
-------
1.97
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
1.76 +
1.56 +
1.35 +
*
1.14 +
.929 +
*
.719 +
.510 + *
o
.301 +
*
.923@-1 +
0 *
* *
0
* * +
o 0 +
* *
+
*
+
+
*
+
~
+
o
+
*
+
* 0
.* * * . * * * ~ ~ . . * * * * *. * * *.
.0 o.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.23
D ER 1 regularization with exact data and
Simpson's rule; A = .150, 8 =1.
77
-------
2.20
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 0 +'
*
1.97 +
1.74 +
1.52 +
1.29 +
1.06 +
.834 + *
.607 +
*
.380 +
0
*
0153 +
0
+
*
*
*
*
+
* * +
o 0
+
* *
+
* +
+
*
+
o
*
+
o
o * * 0
.8 * * e * * * * * * * * * * * * e * * e.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.24
DERI regularization with rounded data and
Simpson's rule; A = .075, 5 = 1.
78
-------
1. 97
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
*
1. 76 +
1.55 + 0
*
1.34 +
*
1.13 +
.920 +
*
.710 +
.499 + *
0
.288 +
.
.778@-1 + 0 *
. * * * . * * * * * .
*
*
*
+
o
*
+
+
*
+
+
*
+
*
+
o
+
*
*
o
. * * * *
+
* . * * *.
.0 O.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.25 DERI regularization with rounded data and
Simpson's rule; A = .150, 0 = 1.
79
-------
2.61
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
2.33 +
2.05 +
1.77 +
1.49 +
1.21 +
.928 +
.648 +
.367 +
.863@-1 +
:~*
+
+
+
+
+
+
+
+
*
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.26 DER2 regularization with exact data and
Simpson's rule; A = .010, e = 1.
80
-------
2.00
1. 79
1. 57
1. 36
1.14
.928
.713
.498
.284
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
*
*
+
*
+
...
+
+
$
(j'j
+
+
*
*
+
+
+
+
*
*
+
+
+
*
+
*
o
o
+
+
*
*
.689@-1 +
.* * * $ *
*
*
+
* Ii * * *.
* . * *
(j'j
$ * * (j'j *
.0 O.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.27 DER2 regularization with exact data and
Simpson's rule; A = .1 00, f = 1
81
-------
3.04
. +. . . . . +. . . . . +. . . . . +. . . . . +. . . . . +. . . . . + . . . . . +. . . . . +. . . . . +. . . . . +. . . . . +. . . . . + . . . . . +. ..
+ +
2.70 +
2.3 7 +
2.03 +
1.69 +
1.36 +
1. 02 +
.681 +
.344 +
.715@-2 ~~.
+
+
+
+
+
+
+
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.28 DER2 regularization with rounded data and
Simpson's rule; A = .010, f = 1.
82
-------
2.06
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 0 +
*
1. 84 I-
1.62 +
.
1.39 +
1. 17 + *
.948 +
*
.725 +
.507 + *
o
.279 +
*
.555@-1 + 0 *
. * * * . * * * * * .
*
*
+
*
*
+
«I
+
*
+
+
*
+
*
+
o
+
*
*
+
* «I * * *.
«I * * " *
.0 O.
.+.... .T.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4~80
Figure 7.29 D ER2 regularization with rounded data and
Simpson's rule; A = .1 00, e = 1.
83
-------
2. 91
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
+
2.61 +
2.32 +
2.02 +
1. 73 +
1. 43 +
1. 14 +
.846 +
.552 +
.257 +
+
+
+
+
+
+
+
*
+
* *
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.30 CDR regularization with exact data and
Simpson's rule; A = O.
84
-------
1. 9 7
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
1.77 +
1.56 +
1. 36 +
1.15 +
.950 +
*
.745 +
.540 +
.
.335 +
*
.130 +
*
* *
* * +
+
* *
0 0
+
* * +
+
"
+
+
$
+
*
+
*
o 0 0 0
.$ * * * * * $ * * * * * * $ * * * * * $.
. +. . . . . +. . . . . +. . . . . +. . . . . +. . . . . +. . . . .+. . . . .+. . . . .+. . . . . +. . . . . +. . . . . +. . . . . +. . . . . +. .
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.31 CDR regularization with exact data and
Simpson's rule; A = .025.
85
-------
3.56
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 0 +
+
'"
'"
+
'"
+
'" '"
+
'"
+
'"
+
'"
'"
+
3.20 +
2.84 +
2.48 +
2.11 +
1. 7 5 +
1. 39 +
1. 03 +
.675 +
'"
.315 +
'"
+
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.32 CDR regularization with rounded data and
Simpson's rule; A = 0
86
-------
1. 98
1.78
1.57
1. 36
1.15
.947
.740
.533
.326
.118
.+. ... .+. . . . .+... ..+. . ...+. .. ..+.. . ..+.... .+.. .. .+. . ...+. . . . .+. .. . .+. ....+.. . . .+..
+ $ +
*
*
+
*
*
+
+
+
*
*
o
o
+
+
+
*
*
+
+
+
*
+
+
+
+
~
$
+
+
*
"
+
+
*
*
.* * * . *
* $ * *
$
. * * . *
* $ * * *.
.0 o.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.33 CDR regularization with rounded data and
Simpson's rule; A = .030.
87
-------
1.97
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ ~ +
*
*
1.77 +
1.57 +
1.37 +
1.17 +
.975 +
*
.775 +
.575 +
*
o
.375 +
*
.175 +
*
*
*
+
+
e
e
+
* * +
+
*
+
+
*
0
+
*
+
*
.~ * * . * * e * * e e * * $ * * e * * $.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.34 CD R regularization with exact data and rec-
tangle rule; A = O.
88
-------
1.97
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ . +
*
*
1.77 +
1.57 +
1.36 +
1.16 +
.956 +
.753 +
.549 + 0
*
.345 +-
*
.142 +
*
*
+
*
o
*
o
'-
*
+
*
+
* *
+
o +
*
+
*
+
* 0
o *
.. * * * * *. * * . . * * . * * * * * ..
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.35 CDR regularization with exact data and rec-
tangle rule; A = .025.
89
-------
2.37
2.13 +
1. 89 +
1.65 +
1. 41 +
1.17 +
.930 +
.690 +
.450 +
.210 +
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
+
+
+
+
+
+
+
+
+
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
Figure 7.36 CDR regularization with rounded data and rec-
tangle rule; A = O.
90
-------
1.97
.+.... .+.... .+.... .+.... .+.... .T.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 8 +
*
*
1.77 +
1.57 +
1.37 +
1.16 +
.962 +
.760 +
.557 +
6)
.355 +
*
.152 +
*
*
+
+
*
*
o
o
+
*
*
+
+
*
*
+
o *
.8 * * * * * 8 * * 8 8 * * 8 * * * * * 8.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
-6.00 -4.20 -2.40 -.600 1.20 3.00 4.80
+
8
+
*
+
* 0
Figure 7.37 CDR regularization with rounded data and rec-
tangle rule; A = .030.
91
-------
Section 8.
Application to Ge(Li) Spectrometer Data
The final problems to which the methods developed under this contract
were applied dealt with the attempt to achieve resolution improvement of
actual spectral records recorded by a Nuclear Data ND100 spectrometer during
an analysis of radium samples.
In conducting these computations it was
assumed that the model of (2.1) is valid and that h can be modeled as a
difference kernel, as previously mentioned in the Summary.
The recorded
response of the instrument to a single radium peak was therefore accepted
as comprising samples of the impulse response of the spectrometer; this
record will hereafter be referred to as IPT4".
The data used in the first
series of these calculations was constructed by superposing 5% random
noise on the PT4 data.
(Throughout this discussion the records which are
described all consist of 31 data points.
In this section all curves are
shown in normalized form.)
This record was given the notation IPT4N05";
it is referred to as "spectral data" in the figures in this section.
By
dealing with data generated in this way it was possible to initially avoid
any tendency of the spectrometer to violate the shift-invariant hypothesis
and the resulting difficulty associated with such pathology of the instrument.
The results of the singular value analysis of the spectral data are
shown in Figures 8.1 through 8.7.
In these figures it can be observed that
values of k large enough to yield any significant resolution improvement in
the data yield also accompanying contamination by noise effects which is
unacceptably large for many purposes.
This situation is typical of the
results obtained in all calculations with actual analyzer data using the
singular value analysis approach.
92
-------
Figures 8.8 through 8.12 indicate the type of results obtained with
Euclidean regularization and various values of the smoothing parameter A.
Here some improvements over the results obtained using singular value
analysis are noted.
For example, Figure 8.11 (A = .001) reveals a fairly
significant level of resolution improvement with the degrading effects
1 imited primarily to some "undershooting" of the zero count level in the
immediate vicinity of the peak tip location.
The other figures indicating results obtained using the spectral data
(Figures 8.13 through 8.21) show typical behavior of the differential and
continuous/discrete regularization methods as applied to spectrometer data.
No clear superiority of these methods (over Euclidean regularization) is
apparent in these results and, due to the larger costs incurred with their
use, there seems to be little to recommend them (in this type of problem).
The lower quality of the results observed to this point in this
section (as compared with the previous one) seems to be due to an inherent
instability in the process of attempting to deconvolve the analyzer peak
data so as to produce the (theoretically correct) extremely sharp peak
comprising the presumed incident data signal.
One knows, of course, that
deconvolution is normally an unstable process in any case, but the results
obtained thus far (compared to those produced in Phillipsl problem, for
example) indicate that the problem at hand is especially unstable.
The
following compromise problem formulation was suggested by C. Nelson,
Environmental Analysis Division, Office of Radiation Programs, as previously
mentioned.
With ~he viewpoint that the impulse response k of the analyzer
is of an approximately Gaussian nature, one may select a Gaussian pulse p
of half-width significantly less than that of k.
Then a resolution improve-
93
-------
ment calculation which yields the result p (with low superposed noise) from
the spectral data would be a highly acceptable, though not ideal, achieve-
ment.
Continuing in the above vein, one may write (2.1) as g = k*f where the
star (*) denotes the operation of convolution.
Convolving with p gives
g*p = k*f*p or
k*f = g
P P
(8.1)
where, through the rest of this section, a subscripting with the name of a
function will be used to indicate the result of convolution with that
function.
(Since convolution is commutative there is no need to insist on
a given order in which the operands are listed.)
If one replaces the pro-
blem of solving (2.1) for f with that of solving (8.1) for fp then the
compromise mentioned above has been conceptually effected. From the
theoretical point of view one may regard the compromise problem in the
following manner.
The problem of solving (2.1) for f (as stated) has proven
to be a quite ill-conditioned one so that significant levels of resolution
improvement are being found extremely difficult to obtain (in the absence
of such a priori information as will be mentioned momentarily) without un-
acceptable contamination due to noise effects.
Therefore one has now re-
placed the original problem with that of determining the moment f given
p
by
v
fp(E) = ~ f(x)p(E-x)dx.
o
is that the solution of (8.1)
(8.2)
The tacit hope, of course,
for the moment f
p
defined by (8.2) will prove to be a significantly more stable problem (when
discretized for actual computation, as usual) than that of solving (2.1) for
f, as originally posed.
94
-------
The pulse p chosen for the data being considered in this section is a
Gaussian of standard deviation 4/3.
The convolution required to produce the
data gp for the problem of (8.1) was, of course, implemented in a discrete
fashion. Although it is perhaps possible to obtain some reduction in com-
putational costs on the ADP, Inc. system through the use of the Fourier
approach to convolution outlined in Section 1, the savings would not be
significant enough in the present problem to justify the additional time
required for software development and so the straight-forward approach was
taken. Specifically, a matrix PM with entries given by PM.. = p(E~E.) was
lJ 1 J
generated; then the product of PM and the vector 9 consisting of samples
~
from g yields the sample vector gp appropriate for a discrete implementation
of (8.1). The approximate solution of the resulting discrete version of
(8.1) was then approached using the Euclidean regularization method.
(Other methods of solution were not considered due to cost considerations
and the prior indications that no significant improvement in the results
obtained is likely.)
When the data g consisted of the previously considered spectral data
~
PT4N05 the g resulting was referred to as the "spread spectral data"; the
p
results obtained are indicated in Figures 8.22 through 8.26. A comparison
of these plots with the previously obtained Figures 8.8 through 8.12 (in
which Euclidean regularization was used with the same spectral data) shows
that very significant improvement of the stability of the problem has been
achieved by using the data spreading technique while only very minimal
losses of resolution level have been incurred.
The remaining consideration of interest in this section has to do
with the application of these methods to other actual photopeak data taken
95
-------
from the analyzer (as opposed to the previous spectral data which was
generated by superposing random noise on the impulse response IPT4" of the
ana lyzer) .
This second data, which will be called the "cross spectral II
data, resulted from the recorded response of the spectrometer to a second
radium peak; it will be denoted as IPT3".
For reasons previously indicated,
only the Euclidean regularization method was applied to the approximate
solution of equations (2.1) and (8.1) with the data g = PT3.
(The impulse
response k = PT4 and the choice of the Gaussian pulse p remained as before.)
The results obtained with this cross spectral, and spread cross
spectral, data are shown by a comparison of Figures 8.27 through 8.31 with
Figures 8.32 through 8.36, respectively.
The greatly increased stability
of the latter group of solutions is most vivid and, in particular, Figure
8.34 shows that acceptable quality results are beginning to become avail-
able through the combined use of the regularization and data spreading
techniques, whereas the regularization methods alone (Figure 8.29) do not
yield results of stability sufficient for many purposes.
One further con-
sideration of gamma spectral data resolution improvement was studied, how-
ever, and found to be of great potential.
One knows, of course, that the solutions f of the problems being con-
sidered in this section are necessarily nonnegative ones.
It has become
obvious, however, that the effects of small amounts of noise and large
tendency toward instability combine to give computed results fA having
significantly large negative values.
It is therefore very reasonable to
hope that one might obtain some improvement in the results of solving both
(2.1) and (8.1) by using the constraint
f(x) ~ 0 for all x.
Computation-
ally this was effected by reconsidering the problem of (2.1) in the
96
-------
Euclidean regularization framework where the a priori information of non-
negativity was now imposed through the appropriate constraints.
statement of the problem of (5.1) was modified to:
-+ -+2 -+2
Minimize Ilg - Afll + A211fll subject to f. ~ 0
E E J
In solving (8.3) computationally an approach analogous
Thus the
(8.3)
to that indi-
cated in Section 5 (following (5.1)) was taken; the well-known subroutine
NNLS of Lawson and Hanson [19] was used in conjunction with a IIdriver
program" written* for the EERF computer system.
Again considering both
the spectral data PT4N05 and the cross-spectral data PT3, resolution
improvement calculations were carried out; the results of these computations
are indicated by Figures 8.37 through 8.40.
These final results appear
most promising.
It is very evident that a great increase in the stability
and sharpness of the enhancements has been obtained in these sample calcu-
lations and, in fact, even setting A = 0 (no regularization at all) appears
to be a completely reasonable alternative.
Using the usual FWHM criterion
one finds that in Figure 8.37, for example, resolution improvement by a
factor of approximately five is in evidence.
Due to lack of time there
were no further calculations conducted with this IIERNNII approach but it is
clear that larger scale investigation of this method should be a high
priority item in the second phase of work under this contract effort.
*Juanita Coley of EERF provided a great deal of expert assistance with
programming and systems analysis during the course of this work.
97
-------
. +. . . . . +. . . . . +. . . . . +. . . . ..... . . . . +. . . . . +. . . . . +. . . . . +. . . . . +. . . . . +. . . . . +. . . . . + . . . . . +. .
.987 + 0 0 +
o '* 0
0 '* 0
.887 + 0 0 +
o
"
0
.787 + 0 +
o '* 0
'*
686 + 0 +
o
0
.586 + 0 +
o
0 '*
.485 .. +
o 0
.3R~ f +
o 0
'*
.284 0 +
o '*
o.
. 184 +0 +
*
'*
.835@-1 + +
'* '*
. '* '* '* '* '* '* '* '* '* * * * * * * '* '*.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1. 00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.1
Singular value analysis with spectral data; k = 1.
98
-------
.984
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .T.... .T.... .T...' .+..
+ 0 $ +
:II
5
.854 + 0
:II
.723 + 0
:II
o
.593 +
*
5
.463 +
0
.333 +
*
.0 0
.203 I-
0 *
.726@-1 + * 0
o *
. * * * * * * :II
-.575@-1 +
0
0
-.188 +
0
+
+
+
+
+
5
+
*
0 +
* 0 0 0 0 o.
* * * * 5 * * * ... *.
o 0 +
o 0
0
+
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.2
Singular value analysis with spectral data; k = 6.
99
-------
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ ~ +
.984
*
*
0
.860 + +
o
.736 + * +
*
.612 +
*
.487 +
0
.363 +
*
.239 +
*
.114 +
0
.0 0 *
0 0 * 0
-.985@-2 +* * II * * * * 0
0 0
--134 + 0
0
+
o
*
+
+
*
+
o
0 0
* +
o
0 * 0
0
* * * * * * * * * *+
o 0
0
0 0+
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.3
Singular value analysis with spectral data; k = 11.
100
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ & +
.866 +
.747 +
.628 +
.510 +
.391 +
.272 +
.153 + *
*
*
+
o
* +
*
+
o
* +
*
+
*
* +
o * 0
.340@-1 + 0 0 0 * 0 0
. * * * * *. * *
.0 0 0 0 0
-.849@-1 +
0
o +
*
0 $ 0 +
& e to * * to * * e e.
o 0 0
+
o
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.4
Singular value analysis with spectral data; k = 16.
101
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
*
*
.868 +
.751 +
.635 +
.518 +
.401 +
.284 +
.167 +
o
+
*
+
*
+
* +
*
0
+
*
* +
o
*
.0 0
.498@-1 + *
o *
. * . . * $ * * 0
o 0
-.671@-1 +
o
+
.
0 0
0
* 0 +
o.
o . . $ * * * *. * *.
o 0
0 0 0 +
o
0
o
o
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.5
Singular value analysis with spectral data; k = 21.
102
-------
.984
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
..
.853 +
.722 +
.590 +
.459 +
.328 +
.197 +
..
..
o
..
o
o
..
o 0
.663@-1 + ..
0 ..
." " .. .. .. * ..
0
-.648@-1 +0 0
o 0
0 0
-.196 +
..
+
+
+
..
+
+
..
0 0 +
o 0
0 "
0
0 0 +
..
0 .. .. .. .. .. .. * .. * ".
o.
o +
o 0
0
+
o
0
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.6
Singular value analysis with spectral data; k = 26.
103
-------
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * * 0 +
.975
0 *
.775 + +
* *
.575 + +
*
*
0 O.
.375 + +
$ 0
*
0 0 0
.175 + 0 * +
*
0 0 * 0
0 0 * * 0 e
-.250@-1 +* * * * * * * * * * * * * * * *+
.0 0 0
0 0
0
-.225 + 0 0 0 0 +
o
0
-.425 + 0 +
-.625 + +
-.825 + +
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.7
Singular value analysis with spectral data; k = 31.
104
-------
1.01
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ ~ 0 +
*
*
.802 +
.591 +
.381 +
.170 +
+
*
*
+
*
*
o
o
o.
+
$
*
o
o
o
o
*
-.251
+
o
o
o
+
* 0
*
0 * * * * * * * * * *.
+
o
0
0
0 +
o
0 +
o 0 *
. * * ~ * * e * * 0
-.404@-1 +0
0
0
o
-.462
+
-.672
+
+
-.883
+
+
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.8
Euclidean regularization with spectral data; A = O.
105
-------
.983
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ e +
o
*
*
*
+
*
* +
+
*
*
+
o
* +
o 0 *
0
+
*
0 0 +
e 0 0 0 0
0 * * * * * * * e e e.
o 0
+
o
0
+
.851 +
.718 +
.586 +
.453 +
.321 +
.188
+
o
.560@-1 + 0
0 *
. * * * * * * *
.0 0 0
-.765@-1 + 0
*
o
o
-.209
+
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.9
Euclidean regularization with spectral data;
A = .00001.
106
-------
.986
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ e +
*
*
.871 +
.756 +
.642 +
.527 +
.412 +
.197 + *
+
*
+
*
o
+
*
+
*
o
+
o
.183 +
*
.681@-1 + 0 *
* 0
.* e * $ e *
.0 0 * 0 0
-.466@-1 + 0 0
+
*
0
+
*
+
$ 0
$ * $ * $ $ $ * $ *.
o 0 O.
o +
000
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.10
Euclidean regularization with spectra] data;
A = .0001.
107
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
.867 +
.749 +
.632 +
.514 +
.396 +
.278 +
.160 + *
*
*
+
o
* +
o
*
+
*
+
*
o
+
o
*
*
+
o
o 0
o
o
+
*
0 0
. +
* * e * * * $ e e *.
o 0 0 0 O.
o
+
o
*
.420@-1 +
.e e
*
* *
o
*
o 0
* *
-.759@-1 +
o 0
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.11
Euclidean regularization with spectral data;
A = .001.
108
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
.869 +
.753 +
.637 +
.521 +
.404 +
*
*
o
o
*
*
+
+
+
o
* 0 +
*
+
*
.288 + +
*
0
.172 + 0 +
*
*
.559@-1 + 0 * 0 0 +
o 0 * * 0 O.
.$ $ * $ * * * * * * I< $ * * $ $ *.
o 0 0 0 0
.603@-1 + 0 0 +
o 0
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.12
Euclidean regularization with spectral data;
A = .010.
109
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 8 +
.866 +
.747 +
'27 +
.508 +
.389 +
.270 +
.150 + *
*
*
o
*
0 *
+
+
+
*
+
*
o
+
*
o
*
+
o
+
*
0 0
* 0 +
o * * 8 * * 8 * 8 8 *.
o 0 0 O.
+
o
.312@-1 +
.8 Gt
8
Gt *
o
8
o
o *
*
*
o
.880@-1 +
o
o
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.13
DERI regularization with spectral data;
A = . 100, [) = 1.
110
-------
.982
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ ~ +
*
.839 + +
*
.696 + * +
.553 + +
*
*
C
.410 + +
o
* 0
.267 + * +
o
* 0
.125 + 0 * +
o 0 0 o.
* * 0
.* * * * * * * 0 0 * * ~ * * ~ * * * *.
-.184@-1 + 0 * +
.0 0 0
0 0 0 0 0 0
-.161 + 0 0 +
o 0
-.304 + +
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.14
DER2 regularization with spectral data;
A = .001, e = 1.
111
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
.865 +
.745 +
.625 +
.505 +
.385 +
.265 +
.145 + *
*
*
o
* 0
*
+
+
+
*
+
*
o
+
o
*
*
+
o
+
*
0
0
$ 0 0 +
* * $ * * $ * * $ $.
o 0 0
+
o
.249@-1 +
.$ $
$
$ *
o
$
o
o *
*
*
o
o
..951@-1 +
o
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.15
DER2 regularization with spectral data;
A = .100, E = 1.
112
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
*
*
0
.869 + 0 +
.752 + * -f
*
.635 + +
o
0
.518 + * +
*
.401 + +
*
.284 + +
*
0
0
.167 + * +
*
0
.504@-1 + 0 0 * * 0 0 +
o * O.
.$. $ $ * * * * * * * e * * $ e *.
o 0 0 0 0
-.665@-1 + 0 +
o 0
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.16
DER2 regularization with spectral data;
'A = .500, f = 1.
113
-------
.986
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .T.... .+.... .+.... .+.... .+.... .+..
+ . +
*
o
.
.870 +
.755 + .
.639 +
.524 + *
o
.408 +
*
.293 +
0
.177 +
*
.620@-1 +0 0 *
0 0 * 0
. * * * * * * *
0
-.534@-1 +
0
0
0 0
+
o
+
*
+
+
o
*
+
+
*
0
+
*
0 0 +
o * 0 0
* * * . * * * * * *.
o
0+
o
0
0 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00
5.50
10.0
14.5
19.0
23.5
28.0
Figure 8.17
DER2 regularization with spectral data;
A = .001, f = 5000.
114
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
.867 +
.749 +
.631 +
.514 +
-396 +
.278 +
.160 + '"
'"
'"
+
o
'" +
o
'"
+
'"
+
'"
o
+
o
'"
'"
+
'"
.418@-1 + 0 0
.$ 8 8 '" '" '" 0 '" '"
0 0 0 '" 0
-.761@-1 +
o
+
'"
0 0
'" +
'" '" $ '" '" '" $ e e.
o 0 0 0
0
+
o
o
o
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.18
CDR regularization with spectral data; A = .001.
115
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
.869 +
.753 +
.637 +
.521 +
.405 +
'"
.288 +
. 1 7 2
+
'"
.561@-1 + 0 '"
o 0 '"
.e e '" '" e '" '" '"
0 0 0
-.601@-1 + 0
o
'"
'"
o
o
'"
'"
'" 0
o
+
+
+
o
+
'"
+
+
'"
0
+
'"
0 0 +
'" 0
'" '" '" '" e '" '" $ e.
o 0 0
0 +
o
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... ..,..... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.19
CD R regularization with spectral data;
A = .010.
116
-------
.986
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... ..,...
+ $ +
*
$
o
.876 +
.766 +
*
0 $
.656 +
.546 +
*
*
.436 + 0
.326 +
*
.216 +
* 0
.105 +
*
0 0 0 *
-.467@-2 +&$ $ * * * $ * 0
+
+
+
+
+
+
*
+
o
*
+
*
0 0 0
0 * * $ * * * $ $ $+
o 0 J 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.20
CDR regularization with spectral data; A = .100.
117
-------
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
.987
* 0
0 *
.884 +
0 0
.781 +
*
*
.678 +
0
.576 +
*
.473 +
0
.370 +
'"
.267 +
0
.164 + *
0
.610@-1 + *
$
./1) & /I) * * * *
+
+
+
o
+
* +
o
+
* +
o
+
* 0
+
*
/I) * * * * '" & & &.
OOO/l)(j 00000
.+.... .+.... .+.... .+.... .+.... .T.... ........ .T.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.21
CDR regularization with spectral data; A = 1.00.
118
-------
.987
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+
-------
.987
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ e +
.882 +
.778 +
.673 +
.568 +
.464 +
*
*
o
* 0
*
+
+
+
+
*
*
+
o
o
.359 +
.254 +
.149 + *
+
*
*
+
o
+
*
0
* +
o 0 0 0
* * e ee e * e e e.
*
.448@-1 +
.(1 e
. e.
e
o e
*
o
o
o
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .T.... .+.... .T.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.23
Euclidean regularization with spread spectral
data; A = .00001.
120
-------
.987
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
.882 +
.777 +
.671 +
.566 +
.461 +
.356 +
.251 +
.146 + *
*
+
o
o +
*
*
+
+
*
o
*
+
o
+
*
*
+
o
f-
o *
* 0 +
o 0 0 0 0
* * $ $ ~ * * * $ $.
*
.407@-1 +
.$ $
o
*
o
* $
$ $
o
~
o
o
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.24
Euclidean regularization with spreac ~pectral
data; A = .0001.
121
-------
.+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . ..+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. .
+ $ +
.986
*
*
.878 + 0
o
.769 +
*
*
.661 +
+
+
+
.552 +
*
o
.444 +
.335 +
*
.227 +
0
*
.118 +
*
0 *
.988@-2 +$$ * * * 0 * 0
0 0 0 * 0 0
0
o +
*
+
+
*
+
u
* +
* 0
$ * $ * * * . $ . *+
o 0 0 0 o.
o
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.25
Euclidean regularization with spread spectral
data; A = .001.
122
-------
.986
.877 +
.768 +
.659 +
.550 +
.441 +
.T. . . . .T. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . - +. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. .
+ Ii +
*
*
o
o
+
+
*
*
+
o 0
+
*
*
+
o +
*
+
* +
o
* 0 0
*. * *. * * & . &+
o 0
0
.332 +
* 0
.223 +
*
.114 +
* 0
0 0 *
.453@-2 +$ * . . * * Ii
0
0 0
000
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.26
Euclidean regularization with spread spectral
data; A = .010.
123
-------
.+. . . . .+. . . . .+. . . . . +. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. .
+ * 0 +
.975
* *
.777 +
*
.579 +
.381 +
0
*
.183 +
* 0
0 * 0 0
0 * @ 0
-.145@-1 +$ * $ * * 0 *
0 0 0
0
o
-.212 +
-.410 +
-.608 +
-.806 +
*
+
+
*
0
+
o
*
0 +
o *
0 0
@ * * *
* * * * * * *+
o
0
0 0 +
o
0+
+
+
o
. +. . . . . +. . . . . +. . . . . +. . . . . +. . . . .+. . . . . +. . . . . +. . . . . +. . . . . +. . . . .+. . . . .+. . . . .+. . . . . +. .
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.27
Euclidean regularization with cross-spectral
data; A = O.
124
-------
.979
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 0 * +
*
*
.811 + 0
*
*
.643 +
0
.476 + *
0 0
0
.308 + 0
*
0
.140 + *
0
* *
.* $ * * * * *
-.277@-1 +
.0
0
-.195 +
0
0
-.363 +
-.531 + 0
0
+
+
*
+
o
+
*
0 0
* 0 +
o
* 0 0
* * * * * * * * * *.
o 0+
o 0
0 +
o
0 +
+
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.28
Euclidean regularization with cross-spectral
data; A = .00001.
125
-------
.+. . . . .+. . . . .+. . . . .+. . . . . +. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. .
+ $ 0 +
.983
0
* *
.844 +
.706 +
*
.568 +
*
.430 +
.291 + 0
*
.153 + 0
*
0 0 *
$
.147@-1 +8$ * * * * . 0 0
-.124 +
0
0
0
-.262 0
+
*
+
+
*
+
o
* +
o
0
+
*
0 0
* 0 0 0
* * * * * * * $ * *+
o o.
o
+
o
0 +
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.29
Euclidean regularization with cross-spectral
data; A = .0001.
126
-------
.984
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ 6) +
* $
.860 + +
*
.735 + 0 +
*
.610 + +
o
.486 + * +
*
.361 + +
*
*
.236 + 0 +
* 0 0
. III + * +
o GI 0
* * 0 0 0
.0 0 * 6) $ * * * 0
-.132@-1 +* * $ * 0 * * $ * * * *+
o 0 o.
o 0
0 0
- .138 + +
o
o 0
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.30
Euclidean regularization with cross-spectral
data; A = .001.
127
-------
.985
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ $ +
*
fI
.866 +
.747 +
.629 +
.510 +
.391 +
.272 +
+
o
*
*
0
*
* 0
+
+
+
+
*
*
+
o
.153 + +
* *
0 0
0 0 * * 0
.346@-1 + 0 * fI * 0 +
.$ * * * * * * * * * * * $ $ fl.
o 0 0 0 0 0 0
0 0
-.843@-1 + +
o
0 0
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.31
Euclidean regularization with cross-spectral
data; A = .010.
128
-------
.985
.+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . . -'-. .
+ 0 +
.865 +
.745 +
.626 +
.506 +
.386 +
.266 +
* $
+
o *
+
*
+
o
* +
o
+
*
+
o
0
0 0 +
*
* 0
* * 0 0 +
o $ * * * * * * *.
o
+
o o.
*
.146
+
*
*
0
.264@-1 + 0 0 0 $
.$ * * k * * 0 0
o o. 0 0
.934@-1 +
o 0
.+. . . . .+. . . . .+. . . . . +. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. .
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.32
Euclidean regularization with spread cross-
spectral data; A = O.
129
-------
.986
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ * +
.871 +
.757 +
.642 +
.528 +
.414 +
.299 +
.185 +
*
*
0 +
o
* +
*
+
o
+
*
*
+
o
* +
*
o
+
*
o 0
.703@-1 + *
o *
.* * * * * * * 0 0
o 0 0
-.441@-1 + 0 0
J<
+
*
* *
* * *
o
* *
* *.
o o.
+
o
o 0
o
o
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.33
Euclidean regularization with spread cross-
spectral data; A = .00001.
130
-------
.986
.+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. . . . .+. .
+ $ +
*
$
.873 + +
.760 + 0 * +
*
.646 + +
o
.533 + +
*
*
.420 + +
o
.307 + +
*
*
.194 + 0 +
*
*
.809@-1 + * 0 0 +
* 0 0
(it 0 * * $ 0 0
.$ $ $ $ $ $ 0 * * $ * * . *.
-.323@-1 + 0 0+
o 0 0 0
o
. +. . . . . +. . . . .+. . . . . +. . . . . +. . . . . +. . . . . +. . . . .+. . . . . +. . . . .+. . . . . +. . . . . +. . . . .+. . . . . +. .
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.34
Euclidean regularization with spread cross-
spectral data; A = .0001.
131
-------
. +. . . . . +. . . . . +. . . . .+. . . . .+. . . . .+. . . . . +. . . . .+. . . . .+. . . . .+. . . . . +. . . . .+. . . . .+. . . . .+. .
+ ~ +
.986
* 0
*
.874 +
0
.761 + *
*
.649 +
0
.537 +
*
*
.425 + 0
+
+
+
+
+
.313 + +
*
*
0
.200 + +
* 0 *
.880@-1 + * 0 +
o * 0
61 0 * * 0 0
.$ $ * * 61 0 * * * * * $ * 61 *.
-.242@-1 + 0 0 0 0 0+
o 0 0 0 0
o
. +. . . . .+. . . . . +. . . . . +. . . . . +. . . . . +. . . . .+. . . . .+. . . . . +. . . . .+. . . . . -r. . . . . +. . . . . +. . . . . +. .
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.35
Euclidean regularization with spread cross-
spectral data; A = .001.
132
-------
.986
.+. . . . .+. . . . .+. . . . .+. . . . . +. . . . . +. . . . .+. . . . .+. . . . .+. . . . .+. . . . . +. . . . .+. . . . . +. . . . .+. .
+ 0 +
o
*
.877 +
.768 +
.659 +
.549 +
.440 +
.331 +
.222 +
+
o
*
+
*
o
+
o
+
*
+
o
+
*
+
*
+
o
* 0
* 0 0
* * * 0 * * 0 0 0+
o 0 0
*
*
o
.112 + *
*
o
0 * 0 * 0
.302@-2 +$ * * 0 * * 0
0 0
0
000
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.36
Euclidean regularization with spread cross-
spectral data; A = .010.
133
-------
.984
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
.884 +
.784 +
.684 +
.584 +
.484 +
.384 +
.284 +
.183 +
*
.833@-1 +
*
*
*
*
+
+
*
* +
+
*
+
*
+
+
*
+
*
+
*
.S S . s. . . 0--0-- 0--0-0--0 .. . .. . .. . e.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.37
ERNN regularization with spectral data; A = O.
134
-------
.964
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
*
*
.866 +
.767 +
.669 +
.571 +
.473 +
.375 +
*
.277 +
.179 +
*
.815@-1 +
*
*
+
+
*
*
+
+
*
* +
+
+
*
+
'"
+
*
0--0-0--0 .. . .. . .. . e.
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.38
ERNN regularization with spectral data; A =
.001.
135
-------
.957
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
* *
.861 + +
.764 + +
*
.667 + * +
.570 + +
.473 + * +
*
.376 + +
.279 + * +
*
. 182 + +
*
*
.848@-1 + * +
*
* *
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.39
ERNN regularization with cross-spectral data;
A = o.
136
-------
.957
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
+ +
* *
.861 + +
.764 + +
*
.667 + * +
.570 + +
.473 + * +
*
.376 + +
.279 + * +
*
.182 + +
*
*
.848@-1 + * +
*
.+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+.... .+..
1.00 5.50 10.0 14.5 19.0 23.5 28.0
Figure 8.40
ERNN regularization with cross-spectral data;
A = .003.
137
-------
BIBLIOGRAPHY
1.
Baker, C.T., Fox, L., Mayers, D.F. and Wright, K., Numerical Solution
of Fredholm Integral Equations of First Kind, Comput. J.
7 (1964)
141-148.
2.
Bendat, J.S. and Piersol, A.G., Random Data:
Procedures, Wiley-Interscience, 1971.
Analysis and Measurement
3.
Blinowska, K.J. and Wessner, E.F., A Method of On-Line Spectra Evaluation
by Means of a Small Computer Employing Fourier Transforms, Nucl.
Inst. and Meth.
118(1974) 597-604.
4.
Bloomfield, P., The Fourier Analysis of rime Series, John Wiley and
Sons, 1977.
5.
Brillinger, D.R., Time Serles:
Data Analysis and Theory, Holt,
Rinehart and Winston, 1975.
6.
Chambless, D.A. and Broadway, J.A., Digital Filtering of Speckle-
photography Data, Experimental Mechanics 19 (1979)
286-289.
7.
Clark, E.L., Jr. and Croll, R.H., Jr., Applications of Digital
Spectral Analysis and Filtering to Aerodynamic Testing, Technical
8.
Report SLA-73-7048A, Sandia Laboratories, Albuquerque, NM.
Cooley, T.W. and Tukey. J.W., An Algorithm for Machine Calculation of
Complex Fourier Series, Math. Compo 19, (1965) 297-301.
9.
Cullum, Jane, Numerical Differentiation and Regularization, SIAM J.
Numer. Anal. 8 (1971) 254-265.
10.
Dines, K.A. and Kak, A.D., Constrained Least Squares Filtering, IEEE
Trans. Acoust., Speech, Signal Processing ASSP-25 (1977) 346-350.
138
-------
11.
Ekstrom, M.P., A Numerical Algorithm for Identifying Spread Functions
of Shift Invariant Imaging Systems, IEEE Transactions on Computers
C-22 (1973) 322-32b.
12.
Golub, G.H. and Reinsch, C., "Singular Value Decomposition and Least
Squares Solutions" in J.H. Wilkinson and C. Reinsch (editors)
Handbook for Automatic Computation, Vol. II:
"Linear Algebra",
Springer, 1971.
13.
Hanson, R.J., A eJumerical Method for Solving Fredholm Integral Equations
of the First Kind Using Singular Values, SIAM J. Numer. Anal. 8
(1971) 616-622.
14.
Hilgers, J.W., On the Equivalence of Regularization and Certain
Reproducing Kernel Hilbert Spaces Approaches for Solving First
Kind Problems, SIAM J. Numer. Anal. 13 (1976) 172-184.
15.
Inouye, T., The Super Resolution of Gamm&-Ray Spectrum, Nucl. Inst.
and Meth. 30 (1964) 224-228.
16.
Inouye, I., Harper, T. and Rasmussen, N.C., Application of Fourier
Transforms to the Analysis of Spectral Data, Nucl. Inst. and
Meth. 67 (1969) 125-132.
17.
Koopmans, L.H., The Spectral Analysis of Time Series, Academic Press,
1974.
18.
19.
Lanczos, C., Applied Analysis, Prentice Hall, 1961.
Lawson, C.L. and Hanson, R.J., Solving Least Squares Problems, Prentice
Hall, 1974.
20.
Perry. W.L., Approximate Solution of Inverse Problems with Piecewise
Continuous Solutions, Radio Science 12 (1977) 637-642.
139
-------
21.
Peterson, G.E. et. al., Singular Value Decomposition and Boron NMR
Spectra in Glass, Journal of Non-Crystalline Solids 23 (1977)
243-259.
22.
Phillips, David L., A Technique for the Numerical Solution of Certain
Integral Equations of the First Kind, J. Assoc. Comput. Mach 9
(1962) 84-87.
23.
Rosenfeld, A., Picture Processing by Computer, Academic Press, 1969.
Rosenfeld, A. and Kak, A.C., Dig1tal Picture Processing, Academic
24.
Press, 1976.
25.
Schaeffel, J., Mullinix, B., Ranson, W. and Swinson, W.F., "Computer
Aided Optical Nondestructive Flaw Detection System for Composite
Materials", U.S. Army Missile Research and Development Command,
Redstone Arsenal. AL, Technical Report T-78-5.
26.
Singleton, R.C., An Algol'ithm for Computing the Mixed Radix Fast
Fourier Transform, IEEE Trans. Audio Electroacoust. AU-17 (1969)
93-103.
27.
Tihonov, A.N.. Regularization of Incorrectly Posed Problems, Soviet
Math. 4 (1963) 1624-1627.
28.
Tihonov, A.N., Solution of Incorrectly Formulated Problems and the
Regularization Method, Soviet Math 4 (1963) 1035-1038.
29.
Tihonov. A.N. and Glasko, V.B., An Approximate Solution of Fredholm
Integral Equations of the First Kind, USSR Compo Math, and Math.
Phys. 4 (1964) 236-247.
30.
Twomey. S., On the Numerical Solution of the Fredholm Integral Equations
of the First Kind by the Inversion of the Linear System Produced by
Quadrature, J. Assoc. Comput. Mach 10 (1963) 97-101.
140
-------
31.
Varah, J.M., A Practical Examination of Some Numerical Methods for
Linear Discrete Ill-Posed Problems, SIAM Review 21 (1979) 100-111.
32.
Wahba, G., Practical Approximate Solutions to Linear Operator
Equations When the Data are Noisy, SIAM J. Numer. Anal. 14 (1977)
651-667.
33.
vJahba, G., "Optimal Smoothing of Density Estimates" in Classification
and Clustering, Academic Press, 1977.
34.
Wahba, G. and Wole, S., A Completely Automatic French Curve:
Fitting
Spline Functions by Cross Validation, Communications in
Statistics 4 (1975) 1-17.
35.
Yule, H.P., Mathematical Smoothing of Gamma Ray Spectra, Nucl, Inst.
and Meth. 54 (1967) 61-65.
141
*U.S.GOVERNMENT PRINTING OFFlCE,1980-647-854/5062.Region 4
-------
EJED EPA 520/5-80-003
Chambless, Donald A.
Radiological data analysis
in the time and. . .
- Due Name and Phone# Mcode
EJED EPA 520/5-80-003
Chambless, Donald A. .
Radiological data analysls
in the time and.. .
Due Name and phone# Meode
,
~_.
'---_.~
1==1--
1__- l
I
U.S. ENVIRONMENTAL PROTECTION AGENCY
Office of Prevention, Pesticides & Toxic Substances
OPPTS Chemical Library (7407)
401 M Street SW
Washington DC 20460
(202) 260-3944
-------
United States
Environmental Protection
Agency
Office of Radiation Programs
E8stern Environ menta I
Radiation Facility
P.O. Box 3009
Montgomery AL 36193
Postage and
Fees Paid
Environmental
Protection
Agency
EPA-G35
~.
-
U.S..AIL
-.
Official Business
Penalty for Private Use $300
Third Class
------- |