A straightforward deconvolution method for use in small computers

A straightforward deconvolution method for use in small computers

NUCLEAR INSTRUMENTS AND METHODS 163 (1979) 5 1 9 - 5 2 2 , (~) N O R T H - H O L L A N D PUBLISHING CO A STRAIGHTFORWARD DECONVOLUTION METHOD FOR USE...

250KB Sizes 1 Downloads 110 Views

NUCLEAR INSTRUMENTS AND METHODS 163 (1979) 5 1 9 - 5 2 2 , (~) N O R T H - H O L L A N D PUBLISHING CO

A STRAIGHTFORWARD DECONVOLUTION METHOD FOR USE IN SMAIJ~ COMPUTERS E SJONTOFT

H C Orsted lnsntut, Copenhagen, Denmark Recewed 6 July 1978 and in rewsed form 25 January 1979 The subject of this paper is the solution of the convolution (or folding) integral, which relates an observed energy spectrum to the true (or onglnal) energy spectrum The techmque presented ns neither of ~terat~ve nature nor does ~t need Fourier transforms of the observed spectrum It works dnrectly on the observed data The method therefore should be statable to small computers

1. Introduction It 1s well known that repeated measurements of a fixed quantity results m a d~strlbuUon of measurements due to the inexactness of the measunng process Often the resulting d~stributlon is the "normal" or the "Gausslan" distribution. When a great number of different quantities are to be measured in the same experiment the resultlng distributions may overlap and it IS often ~mposslble to recognize the locations of the original quantities Many deconvolutlon procedures L) have been made in order to extract ongmal data from measured overlapping distributions Most of them use a great amount of computauonal effort because they work in Fourier space or are lteratlve Consequently they are not very suitable for small computers. This paper presents a method which works in real space and makes a very modest demand on computing capacity and storage The method is based on the folding theorem of Fourier transforms and demands that the Maclaunn series of the reciprocal function of the Fourier transformed convolution function can be calculated 2. Fourier transforms

Before giving the proof it may be useful to state the following concepts from the Fourier transform theory2) ' If f(x) is a function of the variable x, its Fourier transform- if it exists- is defined as the function F(~) = (1/2n) ~

exp(--l~X) f(x)dx

f(x) is deduced by the inverse Fourier transformation f ( x ) = (1/2n) ~

exp0~x) F(c0dc~ -oe

By definition, the folding of two functions, if it exists, is the expression f *h =

f(x-t)

h(t)dt

-a-j

The interest in the folding concept comes from the following theorem:

Theorem: Provided that it exists, the folding of two functions f(x) and h(x) has its Fourier transform (2n) t F(o~)H(oO, an expression-m which F(o0 and H(~z) are the Fourier transforms of f(x) and h(x) respectively; this transformation is reciprocal In order to clarify this the following Fourier transformations are stated without proof y = f f x ) ~ Y = f(cQ,

y = f ' ( x ) ~ Y = I~F(~), 1 Y - tT x/(-~-~) exp ~ - 2 - a J ~

~ , = {1/02a,

ff

Ixl
"

ff

Ixl>a

Y =

=:>y=

exp

/,

(~.n) ~ sm~a

~a

3. A theory of the deconvolution procedure Suppose we have measured energy spectrum g(x) where x represents e.g. the channel number and g(x) the number of counts in the channel x, and suppose we know that this spectrum is the result of

520

E SJONTOFT

an original but unknown spectrum f ( x ) convoluted with the distribution function h(x). The connection between g, f and h is supposed to be

g(x) =

f

4. Two examples

Example 1 : We suppose the original energy spectrum is convoluted with a Gaussian &stribution, that is"

cO

h(u) f ( x - u ) d u

1 exp(_xZ/2a2) h(u) = tr x/(2~)

= h* f

--cO

If we name the Fourier transforms of g, h and f , G(a0, H(o0 and F(o0 respectively we have.

n(ct) = [~-~-u]~ exp(-ct2a2/2)l,

G(ct) = (2~r)½H(ct) F(ct)

and

If H,(~) = 1/H(oO can be expanded in a Maclaurin series

n I (ct) = (2 rt)½exp (ct2e2/2)

Hl(ct) = HI (0) + HI (0) x

From the series expansion of the exponential funcUon we easily find

x ct+H "(o) 1 ~ctz +

+ H(~") (o)



H,(0) = (2~z)~, gl

+

H~3)(0) = 0,

we have F(ct) = G(ct) H, (ct)" (1/2u~ ½,

H';(0) = (2n)icr 2,

H4(0) = (2u)~3tr 4,

,

and from this we get:

and thus because H,(0), HI(0)

f ( x ) = (1/2zr) ½

H',(0) = 0,

,

are constants"

exp(lctX) f(ct)

a 2 ....

a 4-

as

f ( x ) = g ( x ) - - - f g i x ) + --~ g ( a ) ( x ) - -75-.g ( 6 ) ( x ) + ,q,5

--cO

= (1/2x) ½((1/2r0 ½

exp(lctX) Hi(0 ) G(ct) + --cO

+ (1/2r0 ½

exp(lctX) HI(0 ) ctG(ct) + .. ) --cO

= (1/2~) ½[HI(0) O(X) + HI(O) ( - 1 ) g (x) +

H';

+ ~

(0) ( - 1) g"(x) + ] •

Example 2: If we suppose that our spectrum is folded with a box distribution, l e 11/2a, If Ixl a , we have. H(ct) =(~-r~) ~slncta'cta

H i ( a ) = (at0 ¢ aa As f ( x ) is a real function we must have the con&Uon that the s u m of the imaginary parts of the series (ff it is convergent) must be equal to zero If h(x) is an even function we have

f(x)

=

Ha(O) g(x)

-

s i n act

We find in tables 3) that ax a2x 2 7a4x 4 sina-----~ = l + T + - f f - - f f i - +

H'; (0) 2-----T-g"(x) +

+ ~H(14-)(0) 0(4-)(x) -

] (1/2rt)~

+

31a6x 6 3.7-----~+

[a~x~<~].

That is: a2

. .

H,(0) = (2~) ¢, H',(0) = 0, H'~(0) = 3--~ The factor (1/2n) t might have been o m m e d if we had used a nonsymmetric defimtion of the Fourierand the inverse Fourier transform.

7a 4

/4(? ) = 0

Hi')(0)

= 3-5,'

STRAIGHTFORWARD

DECONVOLUTION

And we get a 2 ,, 7a # (4) f(x) = g(x)--~-g (x)+~g (x)31a 6 3- :/i g(6)(x) + It should be noted that the convolution function in both examples was chosen to be area preserving, le" f oo h(x) dx = 1. -oo

5. Numerical example A test spectrum consisting of two Dirac-delta functions was used as an original energy spectrum:

f(x) = 5 0 6 ( x - 1 0 ) + 3 0 6 ( x - 1 4 ) . This spectrum was convoluted with a Gausslan distribution where a = 2. The convoluted spectrum was calculated for 20 equidistant values of x [the histogram (a) in fig 11 The curves (b), (c) and (d) are the deconvoluted energy spectra including respectively the 2nd, 4th and 6th derivative of g(x) approximated by the simple difference formula for second order denvatlves" g"O)=g(j--1)--2g(j)+ + g O + 1), etc T h s means that when, e g , the 4th derivative is included we get the formula (~r= 2)"

f(j) = 2 g ( j - 2 ) -

10g(j-1) + 17O(j)-

- 10#O+l ) + 2g(j+2). I

35

d /

I

I

--

/lc \\ I

30

Curves (b), (c) and (d) are drawn as continuous curves rather than histograms. As might be expected the use of the 6th derivative easily gives oscillations In the solution when only 20x-values are used. We See from fig. lc that a fairly good localisat~on of the two delta functions and a good measure of their respective strength can be obtained when we include the 4th derivative of g(x). As can be seen from the finite difference formula above the signal/noise (S/N) ratio is reduced by a factor 22.3 when the 4th derivative is included using the s~mple difference formula (if the 6th derivative was included, the S/N reduction would be 62 5) It is ewdent that if this deconvolution method is to be used for experimental data with noise present, one has to be very careful when choosing a differentlatlng routine. The differentiation formula used for the spectrum below was made by fitting locally with a third degree polynomml. That is, for each point in the spectrum a third degree Gram polynomlaP) was calculated, w h i c h - i n a least-squares s e n s e - passed between the point and six neighbour points The second order derivative was then calculated from the coefficient of the polynomlum. The fourth- and sixth order derivatives were found in a slmdar way using the second- and fourth order derivatives, respectively Applying orthogonal Gram polynomials for calculating derivatives seems comphcated but is when using the three term recurrence formula (ref 5, p. 240) computationally fast and easy in a small computer. The advantage of the Gram polynomials is the ease of using a weight function for each point in the spectrum. In the example below is seen a spectrum from 137Cs y-rays incident on a 2 5 4 × 7 6 2 × 9 1 5 m m 3 •I06

:/A

25

521

METHOD

35-3O

20

25 15

2

o

10

i 's

5

10 0

5

-5 0

I

5

I

10

I

15

I

20

-

25

Fig ] (a) The convoluted spectrum, (b) deconvoluted spectrum

wtth 2 terms, (c) with 3 terms and (d) with 4 terms

250

, 1 , I , I , I L I , I , I AI LI , I 300

350

400

450

500

ENERGY

Fag 2a Compton spectrum

550

KeV

600

650

700

750

522

E

SJONTOFT

.10 ~

twe m the deconvolut~ons formula and using Gram polynommls to find denvatwes gwes a S/N reducUon of 90 m the energy range near the Compton peaks This factor is to be compared with a S / N reducUon of 1 3 × 105 which would have been the result ff the simple fimte difference formula had been used for the Compton spectrum.

35 30

~ 2s ~

20

6. Conclusion 5 0

250

300

350

400

450

500

550

600

650

700

750

ENERGY KeV

F~g 2b Deconvoluted Compton spectrum

NE 110 plastic scmcdlator from the Ragshospital Whole Body Counter (fig 2a) For such a large scmclllator, multible Compton scattermgs wdl be present, but due to staUsUcal fluctuations m the recording equipment only one broad Compton peak ~s recorded. It is assumed 4) that the folding function ~s Gaussmn and that the standard deviation is pmporUonal to E t, where E is the energy of the channel (channel width is 8.3 keV). In the theory for the deconvolution the standard deviation tr was not allowed to be a function of x, but as far as the vanat~on on x is weak th~s ~s assumed not to ruin the result. Fig. 2b shows the deconvoluted Compton spectrum when a standard deviation tr(E)= 2.5E ½was used One ~s now able to see two Compton peaks. The first peak at 478 keV which corresponds to the maximum energy which can be transferred in a single collision from an incident 662 keV ),-ray and a peak above 525 keV for the double scattered maximum Thin corresponds to values gwen m ref. 4. It seems doubtful whether more mformaUon can be extracted from fig. 2a when the example m fig. ld is considered, which ideally should have been two delta functions Including the 6th denva-

In the method presented above a deconvoluUon techmque is gwen which has the advantage of using only the observed energy spectrum and ~ts denvatwes multiplied by constants. The method seems to be well stated for use m small computers, as the demand of computer program and storage ~s very small. The only problem m using this method ~s connected with the determmauon of the series expansion of the funcuon Hi(a0 when the dlstnbuUon funcUon h(x) is complicated, but this can be done once and for all m a greater computer. The method also gwes an easdy used judgement of how many terms ~t ~s adwsable to use m the deconvoluUon procedure, as one only has to look for wolent oscdlaUons of the higher denvates of g(x) to see where the calculauons should be cut off. The author w~shes to thank the Ragshosp~tal Whole Body Counter Department for providing spectra and S Hoe for his valuable lnformaUon on the experimental results

References l) A F Jones a n d D L Misell, J Phys A Gen Phys 3(1970) 596 2) A Messiah, Quantum mechamcs, vol 1 (North-Holland Publ Co, Amsterdam, 1967) 3) B O Pence and R M Fosters, A short table o/integrals (Gmn and Comp, 1956) 4) G J Hme, Instrumentatton m nuclear medwme, vol 1 (1967) p 589 5) A Ralston, A first course m numerwal analysts (1965)