- Email: [email protected]
The Fourier method in slow neutron time-of-flight spectrometry
N U C L E A R I N S T R U M E N T S A N D M E T H O D S 73
0969) I89-I99; © N O R T H - H O L L A N D
PUBLISHING
CO.
THE F O U R I E R M E T H O D IN S L O W N E U T R O N T I M E - O F - F L I G H T S P E C T R O M E T R Y A. V1RJO
Technical University of Helsinki, Otaniemi, Finland
Received 7 February 1969 The Fourier method, as it is understood here, consists of first measuring the Fourier transform H(f) of the time-of-flight spectrum h(t) in equally spaced frequency points and then forming an estimate for the spectrum itself by Fourier synthesis. The problem of measuring H(f) at discrete frequencies with the help of a periodically modulated incoming beam and a multichannel analyzer is first discussed, and explicite expressions are
derived for a least-squares estimation of H(f) from the analyzer data. The construction of an estimate h(t) for h(t) is then considered, with special emphasis on the "filtering problem", i.e. on getting a resolution function, that has small side lobes but is nevertheless narrow. An expression is also presented for the standard deviation of h(t).
1. Introduction In recent years, a great deal of interest has been focused at two new time-of-flight methods that are less direct than the conventional one, but that can theoretically give a transmission as high as 50%. In the so called correlation method 1-6) the beam intensity is modulated with a periodic binary pseudo-random signal and the time-of-flight spectrum is calculated as the cross-correlation of the modulating signal and the detector output. The Fourier method 7) in turn consists of measuring the Fourier transform of the time-offlight spectrum at several discrete frequencies and forming then the spectrum itself by Fourier synthesis. For this end the beam intensity is modulated with a periodic signal, that has the discrete frequency in question as a basic frequency, and the output of the time-of-flight apparatus is Fourier analysed. The purpose of this work is to investigate the theory of the Fourier method, because this theory has not been presented in the literature. Special attention will be devoted to statistical properties of the estimates that will be formed for the spectrum, and to the filtering of undesired side lobes of the resolution function.
Let us further assume that the system chopper -~ detector can be described by a time-of-flight spectrum h(t), so that if a neutron passes the chopper at the time t', the probability of its being registered by the detector during the time interval (t,t + dt) is h(t-t')dt. The spectrum h(t), when so defined, includes also the flight path between the chopper and the sample. The case o f a polychromatic incoming beam with or without a sample can also be handled with this simple theory, provided that the chopper modulates neutrons of all velocities in the same manner. Starting from the above assumptions and from the assumption of a constant background of intensity u0 it is shown in appendix A that the number of neutrons counted by the detector is a Poisson process with the intensity
z(t)=
H(I) =
Consider an experimental arrangement like that in fig. 1. The neutron beam goes through the chopper via the sample into the detector. Assume that the number of neutrons, that have passed the chopper before the time t, is an inhomogeneous Poisson process 8) with the intensity y(t). This means that the quantities Y1 and }'2 of neutrons that have passed the chopper during two non-overlapping time intervals A a and A2 are independent, Poisson distributed statistical quantities with the mean values
y(t)dt,
( i - 1,2).
(2)
where the background is assumed to be uncorrelated with the input y(t). One way of measuring h(t) is to measure its Fourier transform
2. Measurement of the Fourier transform
E(~';) = f
o Y(t-ta)h(tl)dtl +u°'
o
h(t)e-i"dt,
(3)
point by point with a periodic signal. In eq. (3) we have put ~o = 2 n f for brevity. The inverse transform is
(1)
DETECTOR
CHOPPER
SAMPLE
Fig. 1. Experimental arrangement.
At
189
A. VIRJO
190
h(t) =fl,H(f)e’w’df.
(4)
Let us now assume that the input y(t) is periodic, f being the basic frequency, and let us develop y(t) into a Fourier series y(t) = C,+
f (C,cos(rot)+D,sin(rot}. r=l
The intensity z(t) of the detector periodic, too z(t) = A,+
5
(5)
signal must then be
{A,cos(rot)+B,sin(rwt)}.
(6)
r=l
The Fourier formulas H(rf)
transform
H(f)
can be calculated
from the
(r= 1,2,...),
= (A, - iBr) / (C, - ior),
i H(O) = (&-%)lG.
It is seen that upn = b,, = +_+,“(k) = 0 (p= 1,2 ,...), so that the quantities Apn, B,, and Acp-fjn cannot be measured in this way. In appendix B several procedures are described, by which one can find estimates a,, 6, for the expansion coefficients a,, b, of eq. (9) given a measured set of Z, values. The simple least-squares approach leads to the minimization problem ,ii
(d0+r~1[6ru~(k)+6,u+(k)]-Z~}2=minimum.
If we differentiate with respect to ii,, ii,, 6, and use the orthogonality properties, eqs. (C3)-(C5) of the functions u,(k), u,(k), which are derived in appendix C, we get the equations r ”
(7)
~(-wL”+r+~(-l)
“4-r
= (2/n) C u,(k)Z,,
The following paragraph will be devoted to the measurement of the coefficients A, and B,, when a multichannel analyzer can be used.
E(Z,) = a, + f
CW,(~) +
r=l
h4Wl~
(9)
,... ),
-;( - 1)%-r = (2/n) i: u,(k)% k=l
(r=1,2,...).
i
Let us use an n-channel analyzer that is “rotated” with the basic frequencyfof the chopper signal, eq. (5). The more general situation, where the analyzer is rotated with an integer multiple mf (m= 1,2,...) of the chopper frequency, will be discussed later. The total measuring time will be denoted by T, and it is supposed that T is an integral multiple of l/J Assume that the analyzer is synchronized in such a way that the first channel is open in the time interval (0, 1 /fn), the second channel in the interval (1 /fn, 2/fn) etc. If the number of neutrons that have arrived at the channel k is denoted by Z, (a statistical quantity), a short calculation gives for the mean value of Zk:
k=l
(r=0,1,2
(8)
3. The Fourier analysis of a periodically modulated Poisson process with a multichannel analyzer
(14)
(15)
The summations over v in eq. (15) must be extended over all those values, where the indices are nonnegative. The equations (15) are not independent, it is in fact easy to see that they reduce to exactly n different equations. To prove this, it is only necessary to notice that eqs. (15) do not change, when r is replaced by n + r, and that the equation ii,,, reduces to an identity (for even n). Here the “equation iin,2” means the equation that is derived from the minimum principle (14) by a differentiation with respect to ii,,,. The group of eqs. (15) is thus fully equivalent to the equations ii, )...) ii+1, bi)...) 6,, when n is even, and to the equations ii, ,..., ii,,_+, &i ,..., 6,,_,, when n is odd. As these n equations can be solved for only n of the coefficients ii,, 6,, it is natural to assume in the solution that ii,= g,=O, when r >$z. In this way one gets the following estimates
where (r=O, I,2 ,... ),
a, = (A,T/n)sinc(r/n),
(10)
’
ii0 = (l/n) i
zk,
k=l
(r=1,2,...),
1 b, = (B,T/n)sinc(r/n), u,(k) = cos {(nr/n)(2ku,(k) = sin {(zr/n)(2k sine(x)
(11)
/
(r= 1,2,...),
sin (7rx) = W .
(12)
U,.(k)&,
(1 5 r < in),
k
i
& = (2/d ck u,(k)Zk,
(1 5 r < *a),
I&,, = (l/n) c Q,,(k).&,
(if n is even).
l)}, - l)},
a, = (2/n)z
I
(16)
k
(13)
It can also be directly seen that the estimates
(16) satisfy
SLOW
NEUTRON
TIME-OF-FLIGHT
eqs. (15)approximately, when n>> l, provided that the Fourier coefficients go to zero sufficiently rapidly. The mean values and variances of the corresponding estimates ~4r and/3r are calculated in appendix D. For reference they are written here more explicitly together with the expressions for ~4t,/~ themselves:
Ao = (l/r) ~ z~, k=l
E ( & ) = Ao,
aZ(Ao) = Ao/T,
Ar-
(17)
2 ~ blr(k)Zk' T sinc(r/n) k= 1
E(~r) -= Z r Ap,+ r _ ~, r Ap,_r, p pn+r p pn--r aa(i[r) -
1
{2A o + sinc(2r/n).
r sinc2 (r/n) t pn + 2r
p p n - 2r (18)
2
&-
n
TsinT(r/n)~=1v,(k)Z~,
E(~r)= Z r Bp,+r + Z r Bp,_r, p pn+r p pn-r a2(/~r ) _
1
variance of a linear combination of these estimates, one must start from the more complete expressions (D5, 6). We now turn to the more general case, where the n-channel analyzer is rotated at a harmonic frequency m f ( m = 1,2 .... ). A short calculation shows that the eqs. (9)-(11) are still valid, provided that we replace the quantities at, br, At, Br by arm, brm, Arm, Br,,. It follows that all the formulas and derivations of the present chapter remain true if the indices of all Fourier coefficients and of their estimates are multiplied by m. One conclusion is that, in addition to the previous limitations, it is not possible to get information about A r or B,, if r ~ 0 (rood m). When n is finite and the quantities of primary interest are certain higher harmonics A,, B r with fixed r>ln, then an analysis of eqs. (18) and (19), suitably modified, shows that it is useful to have m = r.
4. Construction of the time-of-flight spectrum by Fourier synthesis We now consider a special type of the modulating signal y(t), that is symmetric both with respect to the origin and to the point t= 1/(2f), so that all sinus terms and all even cosinus terms (except Co) are missing from the expansion, eq. (5). Then the expansion (6) of z(t) necessarily contains only odd harmonic terms (in addition to ,40). An example of such a symmetric signal y(t) is the triangular wave used by WhittemoreT):
y(t)
= Co [ 1 +
2) (cos o,t + cos 3o, t +
{2Ao - sinc (2r/n).
+ 2Ascos 5cot +...)].(21)
TsincZ(r/n)
Because D r = 0 , the estimate ft(rf) of H(rf) can be written
pn +2r AP"+zr - ~p p n ~ 2 r
FI(rf)=(A,-iBr)/C,, (19)
In eqs. (18) and (19) it is supposed that 1 < r < ½ n . The p-summations are to be extended over all those values of p, where the indices keep non-negative. When n is even, the corresponding equations f o r / ~ , are, ( n 1 B,/z Tsinc(½) k=*
Z v./2(k)Zk,
.=o
Ao
az(B"/2) - T sinc 2 (½)"
19i
SPECTROMETRY
(20)
It must be emphasized that the estimates _g~, /3t are correlated, so that if one wants to calculate e.g. the
( r = 1 , 3 , 5 .... ).
(22)
In the estimation of H(0), the background must be taken into account as in eq. (8). In the case n >> 1 the variances of the estimates ~4r,/~, ( r = 1,3,5 .... ) depend only slightly on r. Because the absolute value of Cr decreases rapidly, as r increases, the standard deviation of ft(rf) as calculated from eq. (22) increases rapidly with r. With the incoming signal (21) e.g., the standard deviation of /}(3f) is 9 times that o f / l ( f ) . It follows that the use of the upper harmonics does not greatly increase the information about the system. Assume now that we have made one measurement at each of the frequencies f = A y , 2AI,... , NAy. The analysis of the first harmonic then gives an estimate of H ( f ) at each of these frequencies. In addition, each measurement gives an estimate for H(0). Since all these
A. VIRJO
192
estimates are nearly unbiased, when IZ$1, we may first leave the statistical aspects aside and assume that we know H(f) forf=O, A,, 2A,,.. ., NA,. The question of forming a meaningful estimate for a function h(t), starting from equally spaced values of its Fourier transform H(f), is old and much discussed’, lo). In the following we make especially use of the presentation of Parzen”). We take for the estimate of h(t) h N,~+KJN(~) = A, riNK,(r)ff(rA,)ei”‘r’,
(23)
where o, = 27& and K,,,(r) is a coeficient window which is normalized : KN(0) = 1 and symmetric : KN(r) = KN( - r). We further define a time window k,(z):
We can therefore say equally well that the aliased spectrum h”(t) is seen through the window A,k,(A,t). It is interesting to notice that eq. (27) contains an apparent paradox: the window k,,(t) can be chosen in infinitely many ways without changing the other functions in eq. (26) in other words, the integral equation (26) has infinitely many solutions with respect to kN(~). One way of making the decomposition (25) is simply to cut R,(z) in equal pieces of unit length. In this way one usually gets a discontinuous window. Another method will be presented in the following. We introduce a coefficient window generator K(u) that satisfies the conditions K(0) = 1, K(u) = K( - tf),
k,(z) = r=iNK,(r)e2”i”.
( K(u)=O,when The time window normalized :
n,(r)
is periodic
s
with the period
and symmetric: then be written
= 1
n,(r) = k,( - T). The estimate
h N,+,K&) = Ash(t) * &(A.,4
(23) can
(25)
where the asterisk (*) means the convolution operation. Accordingly we can say that the estimate hN,df,KN(t) shows us the true spectrum h(t) “through the window A,k,(A,t)“. This window is symmetric, periodic with the period l/A, and normalized over the interval (OJ&). Let us write the periodic window position of aperiodic windows :
$
L,(T) =
k,(z-
as a super-
where h(t) is the “aliased h”(t) =
r),
* k,(A,%
(27)
spectrum”
f h(t-r/A,). r=-CC
(30)
Let M be a number depending on N so that: N-C MS N+ 1, which implies that K(r/M) = 0 whenever r=N+ 1, N+2,... . If we now choose the coefficient window generator K(u) in such a way that K(r/M)
(r=O,l,...,
= G(r),
N),
F
G,(z) =
K(r/M)e’“‘*‘.
r=-CC
where kN(r) is symmetric and normalized over the interval (- co, a). An expansion like (26) is always possible, but it is not unique. The criteria that must be followed in the choice of kN(T) naturally depend on, what the estimate of h(t) is needed for. It is often required that kN(5) has a well localized and narrow main peak and that it is attenuated rapidly as z increases. It follows from eqs. (25) and (26), that = A,&)
k(z) = /~mK(u)e2ZirUdU.
(31)
we can write n,(r)
r=-CC
h N,+Y#)
(29)
We assume that K(u) is continuous and has a bounded variation. The time window generator k(z) is then defined as the inverse Fourier transform
1
o fi,(z)dz
lu[zl.
1;
(28)
(32)
The condition (31) only fixes K(u) at the points u = r/M (r=O, f l,..., +N), otherwise K(u) can be chosen arbitrarily. If the Poisson sum formula”) is applied to the function K(u/M)exp(2Cu), one gets
E,(z) =
A4
t
k[M(z-r)].
I=-CC We have thus constructed the following window satisfying eqs. (26) and (27). k,+(T) = Mk(Mr).
(33) aperiodic
(34)
The method described above gives all those aperiodic windows kN(z), whose Fourier transform is continuous. Because the Poisson sum formula is not valid for discontinuous functions, one gets erroneous results if in the above method K(u) is chosen to be discontinu-
SLOW NEUTRON
TIME-OF-FLIGHT
193
SPECTROMETRY
TABLE 1 Five different coefficient windows KN(r ), and the corresponding periodic and aperiodic time windows/~N(r) and kN(r) that are defined by eqs. (24) and (26).
1 Dirichlet : KN(r) = !r l ,
Ku(r)
kN(T)
kN(T)
Irl<=N
/3u(r) _ sin 2~(Nsin 7rr+ ½)3
sin 2~(N=3+ ½)3
0, Irl>N 2
Modified Dirichlet: 1, I r[ < N
/3"(~) sinZ=Uv
KN(r)= ~ ½, [rl=N
sinZ=N~
tan~3
zc~
0, I r l > g 3
Fejer:
KN(r) = 0, 4 Tukey:
[rI>N
rl( l + c o s ~r) , N+I
KN(r) =
I
10,
[rl
(cf. window 2 above) ~~, ~-, [ +
[r[>N
+~DN+ t r
1-,[
1
] +
1]
sin 2~(N + 1)3 27r3[1 - 4(N + 1)232]
2 ( N + 1)
5 Parzen :
~l-6uZ+6[u[ 3, [u[<½ where
K(u)=~ 20-Iul) 3, ~0,
1=<1u1<1
3 ( N + l ) [sin½n(N+l)3 4 4 k ~T)z]
lul>l
ous*. In the case of a continuous K(u) it is in general recommendable to take M = N + l, because one then gets the narrowest ku(z). The table 1 shows five different coefficient windows I 0) KN(r) and the corresponding periodic and aperiodic windows /~N(T) and kN(z). In the first two cases the Fourier transform of the window kN(z) is discontinuous, but is is easy to verify directly that the expansion (26) holds in these cases. The modification in the modified Dirichlet window corresponds to applying to Fourier * This h a s h a p p e n e d in 10) with the ordinary a n d modified Dirichlet windows.
integration the trapezoidal rule instead of averaging with equal weights. One could therefore think that the modified Dirichlet window gives a "better" kN(T) than the ordinary one. Now the functions kN(r ) corresponding to these two windows are equal in form, but the kN(~) corresponding to the ordinary Dirichlet window is a little n a r r o w e r . In this case it is better to compare the periodic windows/~N(r) that are defined uniquely by KN(r), and it can be seen at once that the /~N(T) corresponding to the second window diminishes more rapidly in absolute value between the main peaks than that corresponding to the case 1, especially when
194
A. V)RJO
TABLE 2 Properties of four different window types. Nr~ is proportional to the half width of the time window. The second column indicates the height of the first side lobe (with sign), when the main lobe is normalized to unity.rms is proportional to the standard deviation of the estimate h(t) with all the measuring times Tr chosen to be equal. is same as rms, but with the optimal measuring times, eq. (42).
Window
Dirichlet Fejer Tukey Parzen
Nr~_
First side lobe
0.301/(1+1/2N) -0.215 0.443/(1+ 1/N) 0.046 0.500/(1+ 1/N) -0.025 0.638/(1+1/N) 0.0021
4K)rms
1 0.577 0.612 0.520
1 0.5 0.5 0.375
N is small. If N ~ 1, the second window has practically no advantage as compared to the first one, and it will be left out of discussion in the sequel. The foregoing was only to demonstrate that, because of the arbitrariness in the definition of kN(z), the comparisons must, in the last hand, be made between/~N(r). Table 2 shows some properties of four different window types. The first column shows Nzi, where r½ is defined by the equation kN(z~)= lkN(0 ). If the quantity fmax=NAf is fixed, Nz~ is proportional to the half width of the time window AIkN(Ayt). The other column gives the heighth of the first side peak ofkN(r) in percent of the main peak. It seems to be a general tendency that when the side lobes are reduced, the window gets broader. The last two columns refer to the statistics of the estimate of h(t), and will be discussed later. Let us go back to the case, where the H(rAI) of eq. (23) are estimated from measurements. We therefore assume that one has measured a time T~ at the frequency rAs ( r = 1..... N) by using a signal like that of eq. (21). From these N measurements we can calculate the e s t i m a t e s filo(rAf) and fil~(rAi) with the help of eq. (22). According to the discussion after eq. (22) we can disregard the upper harmonics. If in eq. (23) we estimate H(rAI) by/~a(rAf) and H(0) by
r=l
where the sum of the weights p~ is unity and ~o is an estimate of the background, we get for h(t) the estimate ,~(t) =
In eq. (35) we have written for brevity A l ( r ) = AI(rAI) etc., so that the number in parentheses is an index of the measurement. In order to simplify the analysis of the estimate (35) we shall consider only the limiting case n ~ oo, which implies: a. The estimates Al(r) etc. are unbiased, to that the mean value of h(t) can be calculated from any of eqs. (25) and (27); b. Assuming that y(t) lacks all even harmonics we can conclude from eqs. (D9) and (D10) that the only mutual covariances different from zero between the estimates ~o(r), .,~l(r) and/31(r) are
[ 02[io(r)] = An~T,, <' 02[~1(P)]
= o2[nl(F)]
[
covEAo(r),,4,(r)]
=
cov[Ao(r),/~t(r)]
=
2Ao/T,,
( r = l , 2 ..... N),
A,/T,, B,/T,.
(36)
We further assume that the estimate t~o of the background is measured in a separate measurement taking the time fiT, where N
T=ET, r=l
is the total measuring time. Then a2(ho)= Uo/(flT). With the above assumptions the variance of ~(t) can be shown to be
a2[~(t)] =
~A°A~ 2
flAoTU°r=lLZrrpr2 JI-
+4(C~°l° ) L
+
~Ku(r)[Al(r)c°scnr t+
r=l
L
An
(35)
I K2(r)} "
(37) It seems natural to choosep, = T/T, because this choice minimizes the other term in the expression (37). If we furthermore write the third term in integral form, eq. (37) gets the form P
KH(O)+p]
N
+ (2Ay/C1) 2 KN(r)EXl(r)c°s °)d + Bl(r) sin e)d].
Ao
+Bl(r) sinoo,t] + 8 (Co) 2 L
AoA [ 1+
(Af/Co) [- ~o+ r L= l P/{o(r)l +
=
H(0) + p
,=,
SLOW NEUTRON
TIME-OF-FLIGHT
where p = uo/Co is the "relative background". If N>> 1 and if KN(r) and T, have been chosen somewhat intelligently, then the last term in the square brackets is of the order of N z. It is easy to see that the previous term is at most of the order of N, and it can therefore be disregarded. The same is also true for the second term, provided we have fl>>l/N g. So we can write approximatively
where
1 N T
(4o)
It should be noticed that CoT is the total number of neutrons that have passed the chopper during the total measuring time. The parameter S is of the order of unity and depends, except on the window KN(r), also on the measuring times Tr. We shall consider two ways of choosing the Tr: a. The measuring times are chosen equally long: Tr = TIN. One gets S = ( K >rms = r=l
KZu(r)
T, = KN(r
)/[
)] ,
(42)
KN(r).
(43)
Ku(r r
whence S = = ~
and with a fixed window type it is directly proportional tOfmax. With a given window type we accordingly have the "uncertainty relation" t_~a[-h(t)] = const.,
The table 2 shows the values of ( K ) and (K)rms for different windows. They have been calculated by replacing the summations in eqs. (41) and (43) by integrals, and hold therefore approximately, when N~I. The properties of the estimate h(t) of eq. (35) (with pr=Tr/T) can be condensed in the following three points: 1. The mean value of h(t) is hN,3e,KN(t),which is the convolution of the spectrum h(t) and the periodic window AyfcN(Ait ), or equally well the convolution of the aliased spectrum and the aperiodic window AfkN(Afl). If the window type is fixed and N>> 1, then the half width of the both time windows are inversely proportional to the highest measuring frequency fmax" 2. The standard deviation of h(t) is given by eq. (39),
(44)
where t~ is the half width of the periodic or aperiodic time windows. A "second uncertainty relation" can be seen from table 2: if the standard deviation of h(t) is reduced by the choice of the window type, the main peak gets broader. This relation is, however, only qualitative. 3. The choice of the frequency difference AI determines the amount of overlapping according to eq. (28). It should be mentioned that the window formalism can also be developed in another, fully equivalent mannerg). There one first forms a raw estimate ho(t) of the spectrum by using the Dirichlet w~ndow, and then applies a smoothing operation - a "filter" - on ho(t). For the Tukey window e.g. this smoothing operation takes a very simple form:
= ¼ho t + 2(N+I)A
+
+¼ho [t
(41)
b. The measuring times are chosen to minimize the variance of h(t):
195
SPECTROMETRY
hD(0+
1
]-
2 ( u ; 1)A ,
(45)
The coefficient window formalism of this chapter seems, however, to be simpler in general.
5. Conclusions The window formalism of the previous chapter offers a very flexible method for forming different estimates ofh(t), when the values of H(f) are measured in equally spaced frequency points. Table 2 shows e.g. that by using the Parzen window instead of the Dirichlet window the standard deviation of the estimate can be reduced to 37.5% and the first side lobe to 1% of the corresponding Dirichlet values [-optimal measuring times, eq. (42)]. The approach of the present work can also be used to analyze a cross-correlation measurement. The Fourier spectrum of the pseudo-random signal ordinarily used in cross-correlation is composed of discrete, equally spaced lines, whose amplitude is constant for low frequencies. One gets therefore from a single measurement the equispaced values of H(f) that are needed in order to construct h(t) by Fourier synthesis. This approach is not so simple as the ordinary crosscorrelation, but it is at least more flexible in that the window can be chosen very freely. It is of course possible to apply Fourier analysis to other kinds of
196
A. V I R J O
signals as well, e.g. to the conventional, periodically pulsed signal. An investigation along the above lines is in progress with emphasis on the pseudo-random signal.
It is seen that Z(A 1 I t) and Z(A2 [ t) are independent and Poisson distributed, with the mean values 2pl(t ) and
The author wishes to express his gratitude to P e k k a Jauho for useful discussions, and for the Finnish A E C for financial help.
(t,t+dt) is Poisson distributed with the mean y(t)dt. The total n u m b e r Z(AI) of neutrons coming into the
Appendix A THE STATISTICAL NATURE OF THE DETECTOR OUTPUT
Suppose that neutrons are coming with the intensity
y(t) into a system characterized by the time-of-flight spectrum h(t), in other words the n u m b e r of neutrons that arrive before the time t is a Poisson process with the intensity y(t). The background will be left out of discussion for the moment. Consider two "channels" A~ and 12, i.e. two disjoint sets of the time axis, with b o u n d e d measure. We shall prove that the output of the detector is a Poisson process, whose intensity z(t) can be calculated in a natural way from y(t) and h(t). Suppose that a neutron comes into the detector at the time t. The contribution of this neutron in the channel A~ is Z(AI] t, 1) ( i = 1,2). The characteristic function of the two-dimensional r a n d o m variable Z(Ai] t, 1) is E(exp [iuZ(A~ ] t,1)+ivZ(A 2 ] t,1)]) =
= px(OeiU+pz(t)e~°+q(t),
•/
h(t'-t)dt',
k q(t)
The n u m b e r of neutrons coming in during the interval
channel A i is a r a n d o m variable that is a sum over the contributions of all the intervals dt. Taking into account the addition property of the Poisson distribution, it can be concluded that Z(A a) and Z(A2) are independent and Poisson distributed with the means (i = 1,2)
E[Z(Ai) ] =
dty(t)pi(t ) = -oo
f
dr'
( i = 1,2),
1-pl(t)-p2(t).
(12)
(14)
The proposition stated in the beginning of this appendix has therefore been proved, and eq. (2) follows immediately from eq. (A4), when a background u o is added.
Appendix B SOME ESTIMATION PROCEDURES
Consider n independent and Poisson distributed r a n d o m variables Zk(k= l .... ,n), whose mean values are functions of some parameters ~,~
E(Zk) = Ek(Yi).
(B1)
The estimates for ?~ are derived by solving the minim u m problem
i (Zk--Ek) 2 = min.
(B2)
k=l
A derivation with respect to ?i gives the " L S equations"
q(t)] m, i (Zk--Ek)@Ek/~Yi) = 0,
in other words Z(A i I t,m) is multinomially distributed. Furthermore, if the n u m b e r of neutrons that come in at the time t is Poisson distributed with the mean 2, the corresponding contribution Z(Ailt) is a twodimensional r a n d o m variable with the characteristic function
[iuZ(A1 It) + ivZ(A 2 t)] } = r=0 ~
or
a. The method of least-squares
I f m neutrons have arrived at the time t and if we denote by Z(A~I t,m) the contribution of these neutrons in the channel A~, then the characteristic function of Z(Ai] t,m) is
E(exp
dty(t)h(t'- t).
d~
In this work it can be assumed that the dependence on 71 is linear. In the following we will briefly describe four different 12) methods to estimate ?v
Ai
[pl(t) e i" + pE(t) e iv +
(A3)
p2(t).
(11)
where pi(t) = (
= exp [2pl(t) (e i " - l ) ] e x p [2p2(t)(e ~v- 1)3.
e-Z[Pl(t)ei"+p2(t)eiV+q(t)]'=
( i = 1,2 .... ).
(B3)
k=l
b. The Z2 minimum method We now start from the m i n i m u m problem
(Zk-- Ek)2/Ek = min, k=l
which leads to the " M X S equations""
k=l
(B4)
SLOW NEUTRON
TIME-OF-FLIGHT
If we make the approximation I ( Z k + Ek)/Ek ~ 1, which is justified at least when the average number of neutrons is large in every channel, we get the equations
(Zk--Ek) ~Ek= O, k=1
Ek
( i = 1 , 2 .... ).
(B6)
197
SPECTROMETRY
for a given harmonic, say A1, grows worse if this number of terms is chosen too large, although the measured values Z k themselves are then represented better by the expansion (9).
Appendix C
~)i
These equations can also be derived from the principle of maximum likelihood to he presented in the next section. They will therefore be called the " M L equations".
ORTHOGONALITY PROPERTIES OF THE FUNCTIONS
u,(k)
AND v r ( k )
Let us define a function e(x): 1, if x is an integer,
c. The method of maximum likelihood Let us define the likelihood function
e(x) =
(C1) 0, else.
L = P ( Z , ..... Z,l~,) as the conditional probability of the measured values of Zk, given ?~. The estimates are calculated by finding the values of 7~ that maximize L, or equivalently, logL. The likelihood function L has the form
By starting from the easily verified formula
exp(2nirk/n) = ne(r/n),
(C2)
k=l
it is readily seen that
C = f i (Ez~/Zk!)e -Ek.
(B)
k=l
=
The M L equations (B6) are derived by putting the derivatives of logL with respect to ?~ equal to zero.
d. The modified Z2 minimum method The method may be defined e.g. by minimizing the weighted sum of squares
~ = o.
Lk=l
In eq. (C3) we have put for brevity
(Zk--Ek)21Zk •
-(x) = ( - 1)%(x).
k=l
This leads to the equations
(Zk--Ek) ~Ek _ 0, k=
1
Z k
( i = 1,2 .... ).
(C3)
(B8)
(C4)
The sum formulas for trigonometric functions, together with eq. (12) and eq. (C3) give the following "orthogonality properties" for u~(k) and vr(k):
~])i
These equations can also be derived by putting E k = Z k in the denominator of the M L equations (B6). The possibility that Zk=O for some k must be taken into account in some way. Two of the above estimation procedures, namely the LS method and the modified MXS method, lead to a set of linear equations with respect to 7~. The LS equations are in addition also linear with respect to Zk, which makes the calculation of the statistical properties of the LS estimates very easy. Another advantage of the LS method is that in many cases the LS estimates do not depend on how many of the 7~ one supposes different from zero. This is the case when Ai and B~ are estimated from eqs. (9)-(13). It is probable that the nonlinear ML and MXS estimates of A~ and B~ have smaller variances than the LS estimates. The nonlinear estimates depend, however, on the number of terms that are taken in the Fourier expansion, eq. (9), and it is possible that the estimate
=
k--e-J
'
k
=
-
,
(c5)
k
Z
= O.
k
All the equations (C5) are valid even if some of the indexes r,s are = 0, provided that we define uo(k)- 1. In a similar way it is easy to derive generalized orthogonality properties of higher order.
Appendix D THE MEAN VALUES AND COVARIANCES OF THE ESTIMATES mi AND B ;
The mean values of the estimates di, Bi can he calculated in a straightforward manner by using eqs. (9)(11), eq. (16), and eq. (C5). The result is
198
A. WRJO E(do) = Ao.
E(A,) =
A~ ~=o
E(~r ) =
=
r+s
- -
n
Bs ,=o
E(B,/z)
8-
sine(r/n) p_
(0
r-}-s
sine(r/n)
T-
- ~-
-7
(0
'
~ ( - 1 ) s sinc(s+½) B(~ +_~)n, ~=o sine({)
(if n is even).
(D1)
When n - , oo, the formulas (D1) give in the limit
~ E(X,) = A,,
( r = 0 , 1 , 2 .... ), (D2)
LE(/],)
B,,
=
(r=1,2 .... ).
The estimates d,,/], are in other words unbiased in the limit n--* o o . We will next calculate the mutual covariances of the estimates d~, ~ , as for example cov(~,,L)
For this end we calculate (1 <
=
~(d,L)- ~(£)~(dD.
(D3)
r,s < ½n)
cov(&G) = (2In) 2 ~2 I2 u,(k)u~(l)cov(Z~,Z,) k=l
= (2/n) 2 ~2
1=1
u,(k)u,(k)E(Zk).
(D4)
k=l
By using eqs. (9)-(11) and a generalization of eq. (C5) one can derive the first of the following equations
sinc(r/n)sinc(s/n) ~- , T . + ~ - - -n 1 ~ Apsinc(p/n) [8_(r-s+p)
c°v('~'Bs) = -T .=,
L
-~-,T.
~- T
(r+s+p)_
'
(DS)
,
where (1
cov(~o,Z) = ~,=z ° ~
-7-
[ _ (:~_9 _ _ (~_+_91 ~ ~ov(~o.~.)= ~1~Z ,.sine(p/n) ~ cov(do,&) = Ao/r,
L
where (1 _<_r<{n). For the estimate
Bn]2 we
give here only its variance
, (D6)
SLOW NEUTRON T I M E - O F - F L I G H T SPECTROMETRY a2(/~./2) = Ao / {T sinc 2 (½) }.
(D7)
W h e n n--r oo, the e x p r e s s i o n s ( D 5 - 6 ) r e d u c e to
[ COV(~r,Z{~) ----T1 (Als-'l + A,+,), i 1
cov(~,,B~) = y [ B , + ~ + s g n ( s - r ) B i , _ , l ],
(D8)
1
cov(/3~,/3s) = ~ ( A I s - r l - - A , + r ) , where O # r # s # O a n d s g n ( x ) = x / l x [ the o t h e r possible cases we have cov(£,£)
for x:/:O. I n
= 1 (2Ao + a2r) '
1
(r,0),
c o v ( ~ r , n r ) = "~ B2r,
(D9)
1 ( 2 A o _ A2,),
~
c o v ( A o , A , ) = Ar/T,
lcov(,~0,/3r) = Br/T,
( r = 0 , 1 , 2 .... ), (D10) ( r = l , 2 .... ).
199
References 1) j. Bj6rkman, Internal Report SSI-110 (AB Atomenergi, Studsvik, 1963). 2) A. I. Mogilner, O. A. Salnikov and L. A. Timohin, Instr. Exp. Techn. 2 (1966) 22. 3) L. Pfi.1, N. Kro6, P. Pellionisz, F. Szlfivik and I. Vizi, Proc. 4 'h IAEA Symp. Neutron inelastic scattering 2 (Copenhagen, 1968) p. 407. 4) F. Gompf, W. Reichardt, W. G1/iser and K. H. Beckurtz, ibid., p. 417. 5) K. Sk61d, Nucl. Instr. and Meth. 63 (1968) 114. 6) A. Virjo, Nucl. Instr. and Meth. 63 (1968) 351. 7) j. F. Colwell, P. H. Miller and W. L. Whittemore, 4 th IAEA Symp. Neutron inelastic scattering 2 (Copenhagen, 1968) p. 429. s) A. Papoulis, Probability, random variables and stochastic processes (McGraw-Hill, New York, 1965) p. 554. 9) R. B. Blackman and J. W. Tukey, The measurement of power spectra (Dover Publ., New York, 1958). 10) E. Parzen, Time series analysis papers (Holden Day, 1967) p. 190. 11) M. J. Lighthill, Fourier analysis and generalized functions (Cambridge U. P., 1958) p. 69. lz) M. G. Kendall and A. Stuart, The advanced theory of statistics 2 (Griffin, London, 1961).
0969) I89-I99; © N O R T H - H O L L A N D
PUBLISHING
CO.
THE F O U R I E R M E T H O D IN S L O W N E U T R O N T I M E - O F - F L I G H T S P E C T R O M E T R Y A. V1RJO
Technical University of Helsinki, Otaniemi, Finland
Received 7 February 1969 The Fourier method, as it is understood here, consists of first measuring the Fourier transform H(f) of the time-of-flight spectrum h(t) in equally spaced frequency points and then forming an estimate for the spectrum itself by Fourier synthesis. The problem of measuring H(f) at discrete frequencies with the help of a periodically modulated incoming beam and a multichannel analyzer is first discussed, and explicite expressions are
derived for a least-squares estimation of H(f) from the analyzer data. The construction of an estimate h(t) for h(t) is then considered, with special emphasis on the "filtering problem", i.e. on getting a resolution function, that has small side lobes but is nevertheless narrow. An expression is also presented for the standard deviation of h(t).
1. Introduction In recent years, a great deal of interest has been focused at two new time-of-flight methods that are less direct than the conventional one, but that can theoretically give a transmission as high as 50%. In the so called correlation method 1-6) the beam intensity is modulated with a periodic binary pseudo-random signal and the time-of-flight spectrum is calculated as the cross-correlation of the modulating signal and the detector output. The Fourier method 7) in turn consists of measuring the Fourier transform of the time-offlight spectrum at several discrete frequencies and forming then the spectrum itself by Fourier synthesis. For this end the beam intensity is modulated with a periodic signal, that has the discrete frequency in question as a basic frequency, and the output of the time-of-flight apparatus is Fourier analysed. The purpose of this work is to investigate the theory of the Fourier method, because this theory has not been presented in the literature. Special attention will be devoted to statistical properties of the estimates that will be formed for the spectrum, and to the filtering of undesired side lobes of the resolution function.
Let us further assume that the system chopper -~ detector can be described by a time-of-flight spectrum h(t), so that if a neutron passes the chopper at the time t', the probability of its being registered by the detector during the time interval (t,t + dt) is h(t-t')dt. The spectrum h(t), when so defined, includes also the flight path between the chopper and the sample. The case o f a polychromatic incoming beam with or without a sample can also be handled with this simple theory, provided that the chopper modulates neutrons of all velocities in the same manner. Starting from the above assumptions and from the assumption of a constant background of intensity u0 it is shown in appendix A that the number of neutrons counted by the detector is a Poisson process with the intensity
z(t)=
H(I) =
Consider an experimental arrangement like that in fig. 1. The neutron beam goes through the chopper via the sample into the detector. Assume that the number of neutrons, that have passed the chopper before the time t, is an inhomogeneous Poisson process 8) with the intensity y(t). This means that the quantities Y1 and }'2 of neutrons that have passed the chopper during two non-overlapping time intervals A a and A2 are independent, Poisson distributed statistical quantities with the mean values
y(t)dt,
( i - 1,2).
(2)
where the background is assumed to be uncorrelated with the input y(t). One way of measuring h(t) is to measure its Fourier transform
2. Measurement of the Fourier transform
E(~';) = f
o Y(t-ta)h(tl)dtl +u°'
o
h(t)e-i"dt,
(3)
point by point with a periodic signal. In eq. (3) we have put ~o = 2 n f for brevity. The inverse transform is
(1)
DETECTOR
CHOPPER
SAMPLE
Fig. 1. Experimental arrangement.
At
189
A. VIRJO
190
h(t) =fl,H(f)e’w’df.
(4)
Let us now assume that the input y(t) is periodic, f being the basic frequency, and let us develop y(t) into a Fourier series y(t) = C,+
f (C,cos(rot)+D,sin(rot}. r=l
The intensity z(t) of the detector periodic, too z(t) = A,+
5
(5)
signal must then be
{A,cos(rot)+B,sin(rwt)}.
(6)
r=l
The Fourier formulas H(rf)
transform
H(f)
can be calculated
from the
(r= 1,2,...),
= (A, - iBr) / (C, - ior),
i H(O) = (&-%)lG.
It is seen that upn = b,, = +_+,“(k) = 0 (p= 1,2 ,...), so that the quantities Apn, B,, and Acp-fjn cannot be measured in this way. In appendix B several procedures are described, by which one can find estimates a,, 6, for the expansion coefficients a,, b, of eq. (9) given a measured set of Z, values. The simple least-squares approach leads to the minimization problem ,ii
(d0+r~1[6ru~(k)+6,u+(k)]-Z~}2=minimum.
If we differentiate with respect to ii,, ii,, 6, and use the orthogonality properties, eqs. (C3)-(C5) of the functions u,(k), u,(k), which are derived in appendix C, we get the equations r ”
(7)
~(-wL”+r+~(-l)
“4-r
= (2/n) C u,(k)Z,,
The following paragraph will be devoted to the measurement of the coefficients A, and B,, when a multichannel analyzer can be used.
E(Z,) = a, + f
CW,(~) +
r=l
h4Wl~
(9)
,... ),
-;( - 1)%-r = (2/n) i: u,(k)% k=l
(r=1,2,...).
i
Let us use an n-channel analyzer that is “rotated” with the basic frequencyfof the chopper signal, eq. (5). The more general situation, where the analyzer is rotated with an integer multiple mf (m= 1,2,...) of the chopper frequency, will be discussed later. The total measuring time will be denoted by T, and it is supposed that T is an integral multiple of l/J Assume that the analyzer is synchronized in such a way that the first channel is open in the time interval (0, 1 /fn), the second channel in the interval (1 /fn, 2/fn) etc. If the number of neutrons that have arrived at the channel k is denoted by Z, (a statistical quantity), a short calculation gives for the mean value of Zk:
k=l
(r=0,1,2
(8)
3. The Fourier analysis of a periodically modulated Poisson process with a multichannel analyzer
(14)
(15)
The summations over v in eq. (15) must be extended over all those values, where the indices are nonnegative. The equations (15) are not independent, it is in fact easy to see that they reduce to exactly n different equations. To prove this, it is only necessary to notice that eqs. (15) do not change, when r is replaced by n + r, and that the equation ii,,, reduces to an identity (for even n). Here the “equation iin,2” means the equation that is derived from the minimum principle (14) by a differentiation with respect to ii,,,. The group of eqs. (15) is thus fully equivalent to the equations ii, )...) ii+1, bi)...) 6,, when n is even, and to the equations ii, ,..., ii,,_+, &i ,..., 6,,_,, when n is odd. As these n equations can be solved for only n of the coefficients ii,, 6,, it is natural to assume in the solution that ii,= g,=O, when r >$z. In this way one gets the following estimates
where (r=O, I,2 ,... ),
a, = (A,T/n)sinc(r/n),
(10)
’
ii0 = (l/n) i
zk,
k=l
(r=1,2,...),
1 b, = (B,T/n)sinc(r/n), u,(k) = cos {(nr/n)(2ku,(k) = sin {(zr/n)(2k sine(x)
(11)
/
(r= 1,2,...),
sin (7rx) = W .
(12)
U,.(k)&,
(1 5 r < in),
k
i
& = (2/d ck u,(k)Zk,
(1 5 r < *a),
I&,, = (l/n) c Q,,(k).&,
(if n is even).
l)}, - l)},
a, = (2/n)z
I
(16)
k
(13)
It can also be directly seen that the estimates
(16) satisfy
SLOW
NEUTRON
TIME-OF-FLIGHT
eqs. (15)approximately, when n>> l, provided that the Fourier coefficients go to zero sufficiently rapidly. The mean values and variances of the corresponding estimates ~4r and/3r are calculated in appendix D. For reference they are written here more explicitly together with the expressions for ~4t,/~ themselves:
Ao = (l/r) ~ z~, k=l
E ( & ) = Ao,
aZ(Ao) = Ao/T,
Ar-
(17)
2 ~ blr(k)Zk' T sinc(r/n) k= 1
E(~r) -= Z r Ap,+ r _ ~, r Ap,_r, p pn+r p pn--r aa(i[r) -
1
{2A o + sinc(2r/n).
r sinc2 (r/n) t pn + 2r
p p n - 2r (18)
2
&-
n
TsinT(r/n)~=1v,(k)Z~,
E(~r)= Z r Bp,+r + Z r Bp,_r, p pn+r p pn-r a2(/~r ) _
1
variance of a linear combination of these estimates, one must start from the more complete expressions (D5, 6). We now turn to the more general case, where the n-channel analyzer is rotated at a harmonic frequency m f ( m = 1,2 .... ). A short calculation shows that the eqs. (9)-(11) are still valid, provided that we replace the quantities at, br, At, Br by arm, brm, Arm, Br,,. It follows that all the formulas and derivations of the present chapter remain true if the indices of all Fourier coefficients and of their estimates are multiplied by m. One conclusion is that, in addition to the previous limitations, it is not possible to get information about A r or B,, if r ~ 0 (rood m). When n is finite and the quantities of primary interest are certain higher harmonics A,, B r with fixed r>ln, then an analysis of eqs. (18) and (19), suitably modified, shows that it is useful to have m = r.
4. Construction of the time-of-flight spectrum by Fourier synthesis We now consider a special type of the modulating signal y(t), that is symmetric both with respect to the origin and to the point t= 1/(2f), so that all sinus terms and all even cosinus terms (except Co) are missing from the expansion, eq. (5). Then the expansion (6) of z(t) necessarily contains only odd harmonic terms (in addition to ,40). An example of such a symmetric signal y(t) is the triangular wave used by WhittemoreT):
y(t)
= Co [ 1 +
2) (cos o,t + cos 3o, t +
{2Ao - sinc (2r/n).
+ 2Ascos 5cot +...)].(21)
TsincZ(r/n)
Because D r = 0 , the estimate ft(rf) of H(rf) can be written
pn +2r AP"+zr - ~p p n ~ 2 r
FI(rf)=(A,-iBr)/C,, (19)
In eqs. (18) and (19) it is supposed that 1 < r < ½ n . The p-summations are to be extended over all those values of p, where the indices keep non-negative. When n is even, the corresponding equations f o r / ~ , are, ( n 1 B,/z Tsinc(½) k=*
Z v./2(k)Zk,
.=o
Ao
az(B"/2) - T sinc 2 (½)"
19i
SPECTROMETRY
(20)
It must be emphasized that the estimates _g~, /3t are correlated, so that if one wants to calculate e.g. the
( r = 1 , 3 , 5 .... ).
(22)
In the estimation of H(0), the background must be taken into account as in eq. (8). In the case n >> 1 the variances of the estimates ~4r,/~, ( r = 1,3,5 .... ) depend only slightly on r. Because the absolute value of Cr decreases rapidly, as r increases, the standard deviation of ft(rf) as calculated from eq. (22) increases rapidly with r. With the incoming signal (21) e.g., the standard deviation of /}(3f) is 9 times that o f / l ( f ) . It follows that the use of the upper harmonics does not greatly increase the information about the system. Assume now that we have made one measurement at each of the frequencies f = A y , 2AI,... , NAy. The analysis of the first harmonic then gives an estimate of H ( f ) at each of these frequencies. In addition, each measurement gives an estimate for H(0). Since all these
A. VIRJO
192
estimates are nearly unbiased, when IZ$1, we may first leave the statistical aspects aside and assume that we know H(f) forf=O, A,, 2A,,.. ., NA,. The question of forming a meaningful estimate for a function h(t), starting from equally spaced values of its Fourier transform H(f), is old and much discussed’, lo). In the following we make especially use of the presentation of Parzen”). We take for the estimate of h(t) h N,~+KJN(~) = A, riNK,(r)ff(rA,)ei”‘r’,
(23)
where o, = 27& and K,,,(r) is a coeficient window which is normalized : KN(0) = 1 and symmetric : KN(r) = KN( - r). We further define a time window k,(z):
We can therefore say equally well that the aliased spectrum h”(t) is seen through the window A,k,(A,t). It is interesting to notice that eq. (27) contains an apparent paradox: the window k,,(t) can be chosen in infinitely many ways without changing the other functions in eq. (26) in other words, the integral equation (26) has infinitely many solutions with respect to kN(~). One way of making the decomposition (25) is simply to cut R,(z) in equal pieces of unit length. In this way one usually gets a discontinuous window. Another method will be presented in the following. We introduce a coefficient window generator K(u) that satisfies the conditions K(0) = 1, K(u) = K( - tf),
k,(z) = r=iNK,(r)e2”i”.
( K(u)=O,when The time window normalized :
n,(r)
is periodic
s
with the period
and symmetric: then be written
= 1
n,(r) = k,( - T). The estimate
h N,+,K&) = Ash(t) * &(A.,4
(23) can
(25)
where the asterisk (*) means the convolution operation. Accordingly we can say that the estimate hN,df,KN(t) shows us the true spectrum h(t) “through the window A,k,(A,t)“. This window is symmetric, periodic with the period l/A, and normalized over the interval (OJ&). Let us write the periodic window position of aperiodic windows :
$
L,(T) =
k,(z-
as a super-
where h(t) is the “aliased h”(t) =
r),
* k,(A,%
(27)
spectrum”
f h(t-r/A,). r=-CC
(30)
Let M be a number depending on N so that: N-C MS N+ 1, which implies that K(r/M) = 0 whenever r=N+ 1, N+2,... . If we now choose the coefficient window generator K(u) in such a way that K(r/M)
(r=O,l,...,
= G(r),
N),
F
G,(z) =
K(r/M)e’“‘*‘.
r=-CC
where kN(r) is symmetric and normalized over the interval (- co, a). An expansion like (26) is always possible, but it is not unique. The criteria that must be followed in the choice of kN(T) naturally depend on, what the estimate of h(t) is needed for. It is often required that kN(5) has a well localized and narrow main peak and that it is attenuated rapidly as z increases. It follows from eqs. (25) and (26), that = A,&)
k(z) = /~mK(u)e2ZirUdU.
(31)
we can write n,(r)
r=-CC
h N,+Y#)
(29)
We assume that K(u) is continuous and has a bounded variation. The time window generator k(z) is then defined as the inverse Fourier transform
1
o fi,(z)dz
lu[zl.
1;
(28)
(32)
The condition (31) only fixes K(u) at the points u = r/M (r=O, f l,..., +N), otherwise K(u) can be chosen arbitrarily. If the Poisson sum formula”) is applied to the function K(u/M)exp(2Cu), one gets
E,(z) =
A4
t
k[M(z-r)].
I=-CC We have thus constructed the following window satisfying eqs. (26) and (27). k,+(T) = Mk(Mr).
(33) aperiodic
(34)
The method described above gives all those aperiodic windows kN(z), whose Fourier transform is continuous. Because the Poisson sum formula is not valid for discontinuous functions, one gets erroneous results if in the above method K(u) is chosen to be discontinu-
SLOW NEUTRON
TIME-OF-FLIGHT
193
SPECTROMETRY
TABLE 1 Five different coefficient windows KN(r ), and the corresponding periodic and aperiodic time windows/~N(r) and kN(r) that are defined by eqs. (24) and (26).
1 Dirichlet : KN(r) = !r l ,
Ku(r)
kN(T)
kN(T)
Irl<=N
/3u(r) _ sin 2~(Nsin 7rr+ ½)3
sin 2~(N=3+ ½)3
0, Irl>N 2
Modified Dirichlet: 1, I r[ < N
/3"(~) sinZ=Uv
KN(r)= ~ ½, [rl=N
sinZ=N~
tan~3
zc~
0, I r l > g 3
Fejer:
KN(r) = 0, 4 Tukey:
[rI>N
rl( l + c o s ~r) , N+I
KN(r) =
I
10,
[rl
(cf. window 2 above) ~~, ~-, [ +
[r[>N
+~DN+ t r
1-,[
1
] +
1]
sin 2~(N + 1)3 27r3[1 - 4(N + 1)232]
2 ( N + 1)
5 Parzen :
~l-6uZ+6[u[ 3, [u[<½ where
K(u)=~ 20-Iul) 3, ~0,
1=<1u1<1
3 ( N + l ) [sin½n(N+l)3 4 4 k ~T)z]
lul>l
ous*. In the case of a continuous K(u) it is in general recommendable to take M = N + l, because one then gets the narrowest ku(z). The table 1 shows five different coefficient windows I 0) KN(r) and the corresponding periodic and aperiodic windows /~N(T) and kN(z). In the first two cases the Fourier transform of the window kN(z) is discontinuous, but is is easy to verify directly that the expansion (26) holds in these cases. The modification in the modified Dirichlet window corresponds to applying to Fourier * This h a s h a p p e n e d in 10) with the ordinary a n d modified Dirichlet windows.
integration the trapezoidal rule instead of averaging with equal weights. One could therefore think that the modified Dirichlet window gives a "better" kN(T) than the ordinary one. Now the functions kN(r ) corresponding to these two windows are equal in form, but the kN(~) corresponding to the ordinary Dirichlet window is a little n a r r o w e r . In this case it is better to compare the periodic windows/~N(r) that are defined uniquely by KN(r), and it can be seen at once that the /~N(T) corresponding to the second window diminishes more rapidly in absolute value between the main peaks than that corresponding to the case 1, especially when
194
A. V)RJO
TABLE 2 Properties of four different window types. Nr~ is proportional to the half width of the time window. The second column indicates the height of the first side lobe (with sign), when the main lobe is normalized to unity.
Window
Dirichlet Fejer Tukey Parzen
Nr~_
First side lobe
0.301/(1+1/2N) -0.215 0.443/(1+ 1/N) 0.046 0.500/(1+ 1/N) -0.025 0.638/(1+1/N) 0.0021
4K)rms
1 0.577 0.612 0.520
1 0.5 0.5 0.375
N is small. If N ~ 1, the second window has practically no advantage as compared to the first one, and it will be left out of discussion in the sequel. The foregoing was only to demonstrate that, because of the arbitrariness in the definition of kN(z), the comparisons must, in the last hand, be made between/~N(r). Table 2 shows some properties of four different window types. The first column shows Nzi, where r½ is defined by the equation kN(z~)= lkN(0 ). If the quantity fmax=NAf is fixed, Nz~ is proportional to the half width of the time window AIkN(Ayt). The other column gives the heighth of the first side peak ofkN(r) in percent of the main peak. It seems to be a general tendency that when the side lobes are reduced, the window gets broader. The last two columns refer to the statistics of the estimate of h(t), and will be discussed later. Let us go back to the case, where the H(rAI) of eq. (23) are estimated from measurements. We therefore assume that one has measured a time T~ at the frequency rAs ( r = 1..... N) by using a signal like that of eq. (21). From these N measurements we can calculate the e s t i m a t e s filo(rAf) and fil~(rAi) with the help of eq. (22). According to the discussion after eq. (22) we can disregard the upper harmonics. If in eq. (23) we estimate H(rAI) by/~a(rAf) and H(0) by
r=l
where the sum of the weights p~ is unity and ~o is an estimate of the background, we get for h(t) the estimate ,~(t) =
In eq. (35) we have written for brevity A l ( r ) = AI(rAI) etc., so that the number in parentheses is an index of the measurement. In order to simplify the analysis of the estimate (35) we shall consider only the limiting case n ~ oo, which implies: a. The estimates Al(r) etc. are unbiased, to that the mean value of h(t) can be calculated from any of eqs. (25) and (27); b. Assuming that y(t) lacks all even harmonics we can conclude from eqs. (D9) and (D10) that the only mutual covariances different from zero between the estimates ~o(r), .,~l(r) and/31(r) are
[ 02[io(r)] = An~T,, <' 02[~1(P)]
= o2[nl(F)]
[
covEAo(r),,4,(r)]
=
cov[Ao(r),/~t(r)]
=
2Ao/T,,
( r = l , 2 ..... N),
A,/T,, B,/T,.
(36)
We further assume that the estimate t~o of the background is measured in a separate measurement taking the time fiT, where N
T=ET, r=l
is the total measuring time. Then a2(ho)= Uo/(flT). With the above assumptions the variance of ~(t) can be shown to be
a2[~(t)] =
~A°A~ 2
flAoTU°r=lLZrrpr2 JI-
+4(C~°l° ) L
+
~Ku(r)[Al(r)c°scnr t+
r=l
L
An
(35)
I K2(r)} "
(37) It seems natural to choosep, = T/T, because this choice minimizes the other term in the expression (37). If we furthermore write the third term in integral form, eq. (37) gets the form P
KH(O)+p]
N
+ (2Ay/C1) 2 KN(r)EXl(r)c°s °)d + Bl(r) sin e)d].
Ao
+Bl(r) sinoo,t] + 8 (Co) 2 L
AoA [ 1+
(Af/Co) [- ~o+ r L= l P/{o(r)l +
=
H(0) + p
,=,
SLOW NEUTRON
TIME-OF-FLIGHT
where p = uo/Co is the "relative background". If N>> 1 and if KN(r) and T, have been chosen somewhat intelligently, then the last term in the square brackets is of the order of N z. It is easy to see that the previous term is at most of the order of N, and it can therefore be disregarded. The same is also true for the second term, provided we have fl>>l/N g. So we can write approximatively
where
1 N T
(4o)
It should be noticed that CoT is the total number of neutrons that have passed the chopper during the total measuring time. The parameter S is of the order of unity and depends, except on the window KN(r), also on the measuring times Tr. We shall consider two ways of choosing the Tr: a. The measuring times are chosen equally long: Tr = TIN. One gets S = ( K >rms = r=l
KZu(r)
T, = KN(r
)/[
)] ,
(42)
KN(r).
(43)
Ku(r r
whence S =
and with a fixed window type it is directly proportional tOfmax. With a given window type we accordingly have the "uncertainty relation" t_~a[-h(t)] = const.,
The table 2 shows the values of ( K ) and (K)rms for different windows. They have been calculated by replacing the summations in eqs. (41) and (43) by integrals, and hold therefore approximately, when N~I. The properties of the estimate h(t) of eq. (35) (with pr=Tr/T) can be condensed in the following three points: 1. The mean value of h(t) is hN,3e,KN(t),which is the convolution of the spectrum h(t) and the periodic window AyfcN(Ait ), or equally well the convolution of the aliased spectrum and the aperiodic window AfkN(Afl). If the window type is fixed and N>> 1, then the half width of the both time windows are inversely proportional to the highest measuring frequency fmax" 2. The standard deviation of h(t) is given by eq. (39),
(44)
where t~ is the half width of the periodic or aperiodic time windows. A "second uncertainty relation" can be seen from table 2: if the standard deviation of h(t) is reduced by the choice of the window type, the main peak gets broader. This relation is, however, only qualitative. 3. The choice of the frequency difference AI determines the amount of overlapping according to eq. (28). It should be mentioned that the window formalism can also be developed in another, fully equivalent mannerg). There one first forms a raw estimate ho(t) of the spectrum by using the Dirichlet w~ndow, and then applies a smoothing operation - a "filter" - on ho(t). For the Tukey window e.g. this smoothing operation takes a very simple form:
= ¼ho t + 2(N+I)A
+
+¼ho [t
(41)
b. The measuring times are chosen to minimize the variance of h(t):
195
SPECTROMETRY
hD(0+
1
]-
2 ( u ; 1)A ,
(45)
The coefficient window formalism of this chapter seems, however, to be simpler in general.
5. Conclusions The window formalism of the previous chapter offers a very flexible method for forming different estimates ofh(t), when the values of H(f) are measured in equally spaced frequency points. Table 2 shows e.g. that by using the Parzen window instead of the Dirichlet window the standard deviation of the estimate can be reduced to 37.5% and the first side lobe to 1% of the corresponding Dirichlet values [-optimal measuring times, eq. (42)]. The approach of the present work can also be used to analyze a cross-correlation measurement. The Fourier spectrum of the pseudo-random signal ordinarily used in cross-correlation is composed of discrete, equally spaced lines, whose amplitude is constant for low frequencies. One gets therefore from a single measurement the equispaced values of H(f) that are needed in order to construct h(t) by Fourier synthesis. This approach is not so simple as the ordinary crosscorrelation, but it is at least more flexible in that the window can be chosen very freely. It is of course possible to apply Fourier analysis to other kinds of
196
A. V I R J O
signals as well, e.g. to the conventional, periodically pulsed signal. An investigation along the above lines is in progress with emphasis on the pseudo-random signal.
It is seen that Z(A 1 I t) and Z(A2 [ t) are independent and Poisson distributed, with the mean values 2pl(t ) and
The author wishes to express his gratitude to P e k k a Jauho for useful discussions, and for the Finnish A E C for financial help.
(t,t+dt) is Poisson distributed with the mean y(t)dt. The total n u m b e r Z(AI) of neutrons coming into the
Appendix A THE STATISTICAL NATURE OF THE DETECTOR OUTPUT
Suppose that neutrons are coming with the intensity
y(t) into a system characterized by the time-of-flight spectrum h(t), in other words the n u m b e r of neutrons that arrive before the time t is a Poisson process with the intensity y(t). The background will be left out of discussion for the moment. Consider two "channels" A~ and 12, i.e. two disjoint sets of the time axis, with b o u n d e d measure. We shall prove that the output of the detector is a Poisson process, whose intensity z(t) can be calculated in a natural way from y(t) and h(t). Suppose that a neutron comes into the detector at the time t. The contribution of this neutron in the channel A~ is Z(AI] t, 1) ( i = 1,2). The characteristic function of the two-dimensional r a n d o m variable Z(Ai] t, 1) is E(exp [iuZ(A~ ] t,1)+ivZ(A 2 ] t,1)]) =
= px(OeiU+pz(t)e~°+q(t),
•/
h(t'-t)dt',
k q(t)
The n u m b e r of neutrons coming in during the interval
channel A i is a r a n d o m variable that is a sum over the contributions of all the intervals dt. Taking into account the addition property of the Poisson distribution, it can be concluded that Z(A a) and Z(A2) are independent and Poisson distributed with the means (i = 1,2)
E[Z(Ai) ] =
dty(t)pi(t ) = -oo
f
dr'
( i = 1,2),
1-pl(t)-p2(t).
(12)
(14)
The proposition stated in the beginning of this appendix has therefore been proved, and eq. (2) follows immediately from eq. (A4), when a background u o is added.
Appendix B SOME ESTIMATION PROCEDURES
Consider n independent and Poisson distributed r a n d o m variables Zk(k= l .... ,n), whose mean values are functions of some parameters ~,~
E(Zk) = Ek(Yi).
(B1)
The estimates for ?~ are derived by solving the minim u m problem
i (Zk--Ek) 2 = min.
(B2)
k=l
A derivation with respect to ?i gives the " L S equations"
q(t)] m, i (Zk--Ek)@Ek/~Yi) = 0,
in other words Z(A i I t,m) is multinomially distributed. Furthermore, if the n u m b e r of neutrons that come in at the time t is Poisson distributed with the mean 2, the corresponding contribution Z(Ailt) is a twodimensional r a n d o m variable with the characteristic function
[iuZ(A1 It) + ivZ(A 2 t)] } = r=0 ~
or
a. The method of least-squares
I f m neutrons have arrived at the time t and if we denote by Z(A~I t,m) the contribution of these neutrons in the channel A~, then the characteristic function of Z(Ai] t,m) is
E(exp
dty(t)h(t'- t).
d~
In this work it can be assumed that the dependence on 71 is linear. In the following we will briefly describe four different 12) methods to estimate ?v
Ai
[pl(t) e i" + pE(t) e iv +
(A3)
p2(t).
(11)
where pi(t) = (
= exp [2pl(t) (e i " - l ) ] e x p [2p2(t)(e ~v- 1)3.
e-Z[Pl(t)ei"+p2(t)eiV+q(t)]'=
( i = 1,2 .... ).
(B3)
k=l
b. The Z2 minimum method We now start from the m i n i m u m problem
(Zk-- Ek)2/Ek = min, k=l
which leads to the " M X S equations""
k=l
(B4)
SLOW NEUTRON
TIME-OF-FLIGHT
If we make the approximation I ( Z k + Ek)/Ek ~ 1, which is justified at least when the average number of neutrons is large in every channel, we get the equations
(Zk--Ek) ~Ek= O, k=1
Ek
( i = 1 , 2 .... ).
(B6)
197
SPECTROMETRY
for a given harmonic, say A1, grows worse if this number of terms is chosen too large, although the measured values Z k themselves are then represented better by the expansion (9).
Appendix C
~)i
These equations can also be derived from the principle of maximum likelihood to he presented in the next section. They will therefore be called the " M L equations".
ORTHOGONALITY PROPERTIES OF THE FUNCTIONS
u,(k)
AND v r ( k )
Let us define a function e(x): 1, if x is an integer,
c. The method of maximum likelihood Let us define the likelihood function
e(x) =
(C1) 0, else.
L = P ( Z , ..... Z,l~,) as the conditional probability of the measured values of Zk, given ?~. The estimates are calculated by finding the values of 7~ that maximize L, or equivalently, logL. The likelihood function L has the form
By starting from the easily verified formula
exp(2nirk/n) = ne(r/n),
(C2)
k=l
it is readily seen that
C = f i (Ez~/Zk!)e -Ek.
(B)
k=l
=
The M L equations (B6) are derived by putting the derivatives of logL with respect to ?~ equal to zero.
d. The modified Z2 minimum method The method may be defined e.g. by minimizing the weighted sum of squares
~ = o.
Lk=l
In eq. (C3) we have put for brevity
(Zk--Ek)21Zk •
-(x) = ( - 1)%(x).
k=l
This leads to the equations
(Zk--Ek) ~Ek _ 0, k=
1
Z k
( i = 1,2 .... ).
(C3)
(B8)
(C4)
The sum formulas for trigonometric functions, together with eq. (12) and eq. (C3) give the following "orthogonality properties" for u~(k) and vr(k):
~])i
These equations can also be derived by putting E k = Z k in the denominator of the M L equations (B6). The possibility that Zk=O for some k must be taken into account in some way. Two of the above estimation procedures, namely the LS method and the modified MXS method, lead to a set of linear equations with respect to 7~. The LS equations are in addition also linear with respect to Zk, which makes the calculation of the statistical properties of the LS estimates very easy. Another advantage of the LS method is that in many cases the LS estimates do not depend on how many of the 7~ one supposes different from zero. This is the case when Ai and B~ are estimated from eqs. (9)-(13). It is probable that the nonlinear ML and MXS estimates of A~ and B~ have smaller variances than the LS estimates. The nonlinear estimates depend, however, on the number of terms that are taken in the Fourier expansion, eq. (9), and it is possible that the estimate
=
k--e-J
'
k
=
-
,
(c5)
k
Z
= O.
k
All the equations (C5) are valid even if some of the indexes r,s are = 0, provided that we define uo(k)- 1. In a similar way it is easy to derive generalized orthogonality properties of higher order.
Appendix D THE MEAN VALUES AND COVARIANCES OF THE ESTIMATES mi AND B ;
The mean values of the estimates di, Bi can he calculated in a straightforward manner by using eqs. (9)(11), eq. (16), and eq. (C5). The result is
198
A. WRJO E(do) = Ao.
E(A,) =
A~ ~=o
E(~r ) =
=
r+s
- -
n
Bs ,=o
E(B,/z)
8-
sine(r/n) p_
(0
r-}-s
sine(r/n)
T-
- ~-
-7
(0
'
~ ( - 1 ) s sinc(s+½) B(~ +_~)n, ~=o sine({)
(if n is even).
(D1)
When n - , oo, the formulas (D1) give in the limit
~ E(X,) = A,,
( r = 0 , 1 , 2 .... ), (D2)
LE(/],)
B,,
=
(r=1,2 .... ).
The estimates d,,/], are in other words unbiased in the limit n--* o o . We will next calculate the mutual covariances of the estimates d~, ~ , as for example cov(~,,L)
For this end we calculate (1 <
=
~(d,L)- ~(£)~(dD.
(D3)
r,s < ½n)
cov(&G) = (2In) 2 ~2 I2 u,(k)u~(l)cov(Z~,Z,) k=l
= (2/n) 2 ~2
1=1
u,(k)u,(k)E(Zk).
(D4)
k=l
By using eqs. (9)-(11) and a generalization of eq. (C5) one can derive the first of the following equations
sinc(r/n)sinc(s/n) ~- , T . + ~ - - -n 1 ~ Apsinc(p/n) [8_(r-s+p)
c°v('~'Bs) = -T .=,
L
-~-,T.
~- T
(r+s+p)_
'
(DS)
,
where (1
cov(~o,Z) = ~,=z ° ~
-7-
[ _ (:~_9 _ _ (~_+_91 ~ ~ov(~o.~.)= ~1~Z ,.sine(p/n) ~ cov(do,&) = Ao/r,
L
where (1 _<_r<{n). For the estimate
Bn]2 we
give here only its variance
, (D6)
SLOW NEUTRON T I M E - O F - F L I G H T SPECTROMETRY a2(/~./2) = Ao / {T sinc 2 (½) }.
(D7)
W h e n n--r oo, the e x p r e s s i o n s ( D 5 - 6 ) r e d u c e to
[ COV(~r,Z{~) ----T1 (Als-'l + A,+,), i 1
cov(~,,B~) = y [ B , + ~ + s g n ( s - r ) B i , _ , l ],
(D8)
1
cov(/3~,/3s) = ~ ( A I s - r l - - A , + r ) , where O # r # s # O a n d s g n ( x ) = x / l x [ the o t h e r possible cases we have cov(£,£)
for x:/:O. I n
= 1 (2Ao + a2r) '
1
(r,0),
c o v ( ~ r , n r ) = "~ B2r,
(D9)
1 ( 2 A o _ A2,),
~
c o v ( A o , A , ) = Ar/T,
lcov(,~0,/3r) = Br/T,
( r = 0 , 1 , 2 .... ), (D10) ( r = l , 2 .... ).
199
References 1) j. Bj6rkman, Internal Report SSI-110 (AB Atomenergi, Studsvik, 1963). 2) A. I. Mogilner, O. A. Salnikov and L. A. Timohin, Instr. Exp. Techn. 2 (1966) 22. 3) L. Pfi.1, N. Kro6, P. Pellionisz, F. Szlfivik and I. Vizi, Proc. 4 'h IAEA Symp. Neutron inelastic scattering 2 (Copenhagen, 1968) p. 407. 4) F. Gompf, W. Reichardt, W. G1/iser and K. H. Beckurtz, ibid., p. 417. 5) K. Sk61d, Nucl. Instr. and Meth. 63 (1968) 114. 6) A. Virjo, Nucl. Instr. and Meth. 63 (1968) 351. 7) j. F. Colwell, P. H. Miller and W. L. Whittemore, 4 th IAEA Symp. Neutron inelastic scattering 2 (Copenhagen, 1968) p. 429. s) A. Papoulis, Probability, random variables and stochastic processes (McGraw-Hill, New York, 1965) p. 554. 9) R. B. Blackman and J. W. Tukey, The measurement of power spectra (Dover Publ., New York, 1958). 10) E. Parzen, Time series analysis papers (Holden Day, 1967) p. 190. 11) M. J. Lighthill, Fourier analysis and generalized functions (Cambridge U. P., 1958) p. 69. lz) M. G. Kendall and A. Stuart, The advanced theory of statistics 2 (Griffin, London, 1961).