A deconvolution method for use in small computers

A deconvolution method for use in small computers

532 Nuclear A deconvolution A. Cabral-Prieto, Instituto and Methods in Physics Research B54 (1991) 532-537 North-Holland method for use in smal...

560KB Sizes 1 Downloads 56 Views

532

Nuclear

A deconvolution A. Cabral-Prieto, Instituto

and Methods

in Physics

Research

B54 (1991) 532-537 North-Holland

method for use in small computers

H. Jim~nez-Do~ngu~z

Nacionol de Inuesti~ac~o~es Nwleares,

Received

Instruments

Apartado

and M. Torres-Valderrama Postal No. 18-1027, Cal. Escandin,

M&co

II801

D.F., Mexico

13 July 1990

An extension of an existing method to defold spectra is presented. The expression to recover any broadening distribution folded with a Lorentzian Line is derived on the basis of the Fourier analysis but its computation is done in real space. Both expressions, the one derived here and the one existing in the literature are applied to broadened simulated and experimental Mlissbauer spectra, where both the Lorentzian signal and the broadening distribution were recovered under certain constraints of the method

1. Introduction In spectroscopies such as Miissbauer, EPR, NMR, etc. Lorentzian spectral lines determined by the Iifetime of the quantum excited states are perturbed or broadened by external effects [1,2]. If, for example, such effects have a random character, then the observed line is the convolution of the Lorentzian signal with a Gaussian dist~bution. The spectroscopist frequently needs to defold his spectra in order to know the Lorentzian and broadening contributions. Many defolding methods have appeared in the literature in order to get either the intelligence signal 13-61, or the broadening function [7-10) from the observed spectra, but most of them require a great amount of computational effort because they either work in Fourier space or they are too lengthy. In this article, an existing method [5,6] based on the Fourier analysis is used to deconvolute Mijssbauer spectra in order to recover the Lorentzian profile whenever the Gaussian cont~bution, which broadens the spectra, is small. This method is extended in the present work by deriving an expression to regain the broadening function, whatever that may be. This expression applies whenever the Lorentzian contribution is small. Examples are given for both simulated and experimental Mijssbauer spectra. The results are compared to those obtained when the simulated and experimental data are handled with a precise deconvoluting method.

where Z(x) represents the observed or convoluted line, and 1(x - y) the intelligence signal, e.g. a Lorentzian line which is perturbed or broadened by a function g(y) that might have a Gaussian shape. Sjontoft [5] has presented a defolding method based on the convolution theorem of Fourier transforms, according to which the folding of two functions, Z(x) and g(x), has its Fourier transform Ill], given by Z(a) = (2T)“‘L(a)G(a),

f2)

where I(a), L(a) and G(a) are the Fourier transforms of the observed Z(x), the Lorentzian Z(X) and the Gaussian g(x) lines, respectively. From expression (2), one can get either the Fourier tr~sform L(a): L(a)

= (27r)-“‘Z(a)G1(a),

with G,(a)

(3) or the Fourier transform G(a): G(a)

= (2-n)-“‘I(a)L,(a),

with L,(a)

depending on which function the experimenter wants to recover from the convoluted spectra. Here G,(a) and L,(a) have continuous derivatives of all orders. The intelligence signal Z(x) may be expressed in terms of the inverse Fourier Transform as: f(x)

The convolution of two functions may be written in a general form as

which after substitution of eq. (3) becomes

Z(x) = ~_$Y)Z(~

Z(x) = (2n)-‘~_WWei”XZ(a)G,(a)

0168-583X/91/$03.50

dy>

0 1991 - Elsevier Science Publishers

(I)

= -$, (4)

2. Theory

-Y)

= &,

= (2n)-“2_/_mme’“xL(a)

B.V. (North-Holland)

da,

da.

(5)

A. Cabral-Prieto

et al. / A deconvolution method for small computers

In Sjontoft’s paper a Gaussian distribution was assumed in one of the examples for the perturbing or broadening function g(x). The reciprocal of its Fourier transform G,(a) was expanded in a Maclaurin series and eq. (5) was worked out and simplified to give [5]:

I(x)

=Z(x) +

-


.‘.)

i.e. the intelligence signal Z(x) is expressed in terms of the observed spectrum Z(x), of its even order derivatives and of u, the standard deviation of the Gaussian distribution g(x). The right-hand side of eq. (6) approximates Z(x) well enough whenever r, 2 To, where r,, r, are the FWHMs of the Lorentzian and Gaussian shapes, respectively. In the present article, the opposite case which might be of practical interest is worked out i.e., the case where the broadening function g(x), whatever that may be, is to be recovered. Following Sjontoft’s procedure [S], the broadening function g(x) can be expressed in terms of its inverse Fourier transform as g(x)

= (~IT-“~/~

eiaXG(a)

do,

-m which after substitution g(x)

of eq. (4) becomes

= (2~)-‘J_“_e’““l(~)L,(n)

da.

(7)

Now, assuming a normalized Lorentzian line as the intelligence signal and expanding Lo in a Maclaurin series, one gets for the broadening distribution g(x): g(x)=

E (-1)” n=O

533

3. Numerical calculations and discussion In this work, the Voigt integral is used to simulate convoluted spectra with higher precision [12-141, and to fit Miissbauer spectra. In simulated data, the extracted spectra from the defolding procedure are to be compared with the area-normalized Lorentzian and Gaussian lines which are used as reference. In experimental data, a smoothing algorithm is used with and without background before applying eq. (9). In both simulated and experimental data, eqs. (6) and (9) are handled according to Verma’s work [6]. There, the best recovered spectrum Z(x) or g(x) is a function of NT, and of m, where NT is the chosen number of terms in the series of eqs. (6) and (9), and m is the separation between the data points to evaluate the corresponding derivatives. These even order derivatives of the convoluted spectra are evaluated according to the general difference formula given by Z(“)(i)

= [ ZCn-‘) (i-m)-2Z(“-2)(i) +ZCnm2)( i + m)] /m2,

where i = 1, 2, 3, ..., N, N is the n th data point and 111the step size to evaluate the even order derivatives Z(“-2)( i T m). On the other hand, eq. (6) will be applied to simulated data only, whereas eq. (9) will be applied to both simulated data and to the stainless steel Mossbauer spectrum. 3.1. Voigt-simulated

spectra

The simulated known expression

line Z(x), is proportional to a well called the Voigt profile [12-141

Z(x)

u),

= =&H(a,

(10)

in which

(8) where p = r,/2, i = m. It can be readily seen that when l(x) is an even function, the imaginary part of eq. (8) is zero. Therefore this expression may be written as:

g(x)

$z’yx)+&z(4yx) ri? - mz(6)(x) + ...

= Z(x) -

Analogously to eq. (6), eq. (9) keeps the same dependence on Z(x) and on its even order derivatives, and this time on r, instead of u. The right-hand side of eq. (9) approximates g(x) well enough whenever r, I r,. As can be seen, expressions (6) and (9) depend on the observed spectrum, on its even order derivatives and on some known real constants so that only the real space is required.

m

H(a,

u) = z

eey2d y

/ -m(u-y)2+az’

(11)

where area-normalized Lorentzian and Gaussian functions are used to get this integral with a = (Z’,/Z”)m, and u = (2(v - v&m)/& which is the dimensionless variable representing the frequency of the Voigtian with respect to the centre v,. Z’o and the standard deviation u of the Gaussian are related by e= r,/2m. The simulated convoluted spectra, with any desired ratio r,/r,, are computed making use of the relationship that exists between the Voigt profile and the complex error function w(z) [13,14]: H( a, u) = Re w(iz), with z = a + iv, and a and v as defined above. Two broadened simulated singlet spectra are chosen:

-2.00

00

a

-0.67 RELATIVE

-0.67 RELATIVE

0.67 VELOCITY

0.67 VELOCITY

=

2.00

NT

2.00

NT = m =

2

2

k-

0.03

q IJ q reference

-2.00

Do

-0.67 RELATIVE

-0.67 RELATIVE

0.67 VELOCITY

0.67 VELOCITY

2.00

NT = 3 In =2

2.00

NT = 2 m = 2

with eq. (6); ~ convoluted Lorentzian-Gaussian line obtained with eq. (10); * profile obtained with eq. (9); convoluted Lorentzian-Gaussian line; 0 q q and (d) the differences between the ordinates of curves •I 0 q and * * are heights for NT = 2, m = 2 is 2.56%.

0.00

0.11

2

i! -

c

x-0.06

Fig. 1. These figures depict the simulated spectra. (a) + * Lorentzian profile obtained area-normalized Lorentzian line. (b) labels are the same as in (a). (c) * reference area-normalized Gaussian line. (d) Labels are the same as in (c). Notice that in (c) amplified. The relative difference between the

0.00

0.04

2 ;L )_

?

p 0.07

0.11

0.00

0.03

P ?

E

go.05

0.08

A. Cabral-Prieto et al. / A deconoolution method for small computers

m = 2. As previously mentioned, deviations in the maximum height and base line begin to appear and neither these nor the FWHM are regained. Also, the characteristic pods in this figure are seen to grow as NT and m increase. In Verma’s paper [6], one can observe the sort of oscillations produced for values as high as NT = 4 and m = 3 in a different spectrum.

one to regain the pure Lorentzian profile and the other one to regain the Gaussian profile. Fig. la shows the deconvolution of the simulated singlet spectrum where the Lorentzian line is recovered. Here, the parameters r,, I-, and velocity per channel are given the values: 0.8, 0.4 mm/s and 0.1 (mm/s)/channel, respectively. When eq. (6) is applied, one can observe from fig. la that all the characteristics of the Lorentzian line are recovered. Reference and recovered Lorentzian lines superimpose each other in this figure. Table 1 shows the values of NT and m for which the height of the recovered Lorentzian line is matched up to the third decimal place with that of the reference line. Beyond these values for NT and m, deviations in the maximum height, FWHM and base line of the curve begin to appear as shown in fig. lb. As can be seen from fig. lb, the recovered height of the Lorentzian shape (continuous-crossed line) is larger than that of the reference line (squared line), whereas the recovered FWHM is smaller than the full width of the reference line. For larger values of NT and m, higher oscillations appear, These cases are not shown. Fig. lc shows the deconvolution of a different simulated singlet Lorentzian-Gaussian convoluted spectrum (continuous line), where the broadening function, i.e., the Gaussian line (continuous crossed line), is recovered. Here, the parameters r,, To and velocity per channel are given the values 0.4, 0.8 mm/s and 0.1 (mm/s)/channel, respectively. In this case, eq. (9) has been applied with NT = 2 and m = 2. It can also be seen from this figure that neither the height nor the FWHM of the reference Gaussian (squared line) are recovered completely. On the other hand, the resulting line shows small pods on the wings. All these features are seen to grow when parameters NT and m are increased. Table 1 shows the values of NT and m for which the regained height approximates the reference Gaussian (squared line) up to the second decimal place. Fig. Id and table 1 shows the data of the same singlet Lorentzian-Gaussian convoluted spectrum of fig. lc, but eq. (9) is calculated for values NT = 3 and

Table 1 Results of the simulated

m

3.2. Experimental

The experimental Mossbauer spectrum of stainless steel is used in this analysis. This spectrum is known to show a broader line width than the natural one which is 0.196 mm/s [15]. Also, the width r, of the broadening function is larger than the Lorentzian width r,, so that by the present defolding method one may be able to regain the broadening function when eq. (9) is used. Eq. (9) is applied to this spectrum in two modes: i) The background is calculated and removed from the raw data (squared line, fig. 2a), then these data are smoothed (continuous line, fig. 2a) by employing an algorithm reported in ref. [16] with a smoothing parameter RHO = 0.00025. The purpose of the smoothing procedure is to reduce oscillations in the defolding operation. ii) The raw data are smoothed (continuous line), and eq. (9) is applied to extract the broadening function (crossed line fig. 2b). The extracted broadening function shows an asymmetrical shape with a distorted oscillatory base line. Fig. 2c shows the deconvolution of the backgroundless stainless steel spectrum where the smoothed data are normalized with respect to the area before applying eq. (9). When this equation is used, the output spectrum (crossed line) shows similar characteristics to those of fig. 2b. However, this time the base line has similarities with figs. lc and Id, that is, a characteristic pod appears on the wings suggesting that a Gaussian line may be involved. These pods are observed only when the background is removed from the raw data. Compare this figure with fig. 2b.

m

Gaussian

h,=0.08Sa

h, = 0.117 a

r, = 0.8 mm/s & = 0.4 mm/s

r, = 0.4 mm/s r, = 0.8 mm/s

2

3

4

5

6

7

0.085 0.106

0.086 -

0.087 -

0.087 _

0.087 _

0.088

’ The heights respectively.

of the recovered

spectra

spectra

Lorentzian

NT: 1 2

convoluted

535

lines are to be compared

1 2

2

3

4

5

6

7

0.090 0.114

0.091 0.121

0.091 -

0.091 -

0.091 _

0.091

with the normalized

heights

h, and h, of the Lorentzian

and Gaussian,

86 CHANNEL

CHANNEL

143 NUMBER

143 NUMBER

2

m

=

= 2

NT

= 0.000250

c 30

143 NUMBER

86 CHANNEL

143 NUMBER

n

86 CHANNEL

200

200

=2

=2

Fig. 2. These figures depict the experimental stainless steel Mossbatter spectrum. (a) 0 0 0 Backgroundless stainless steel Miissbauer spectrum which has been smoothed (smoothed stainless steel spectrum. The background is present; * * * + resulting broadening function when eq. (9) is applied. (c) by using the algorithm of ref. [16]. (b) fitted Gaussian. same backgroundless smoothed spectrum of (a); * * * * output broadening function of eq. (9); (d) * * * * same spectrum as in (c); -

8 d

06

RHO

)

LL

R \

a

a

A. Cabral-Prieto

et al. / A deconvolution

Table 2 Analysis of the stainless steel MBssbauer spectrum Profile

rl_

rG

Lorentzian Voigtian

0.397 0.274 (free) 0.196 (fixed) 0.196 (fixed) 0.186

0.229 (free) 0.293 (free) 0.301 (free) 0.228

Eq. (9) = Eq. (2) [I51

x2

1.91 1.89 1.85

a The parameters NT and m are 2 and 2, respectively. To interpret more clearly the resulting broadening distribution (crossed line, fig. 2c) a fitting procedure is applied to the raw data of stainless steel, employing first a Lorentzian profile as a model, from which we get a value for r,. Secondly, the Voigt profile is used to fit those data leaving both parameters r, and To free in a first run; in a second run, parameter r, is fixed to its natural width, leaving r, free. Finally, the resulting broadening function of fig. 2c is fitted with a Gaussian and a value for r, is obtained. See table 2 for results. When using a Lorentzian line as a model, from the fitting one gets a value for the FWHM of 0.397 mm/s, i.e. a residual width of 0.201 mm/s with respect to the natural line width for the gamma transition of 14.4 keV. The residual width is known to be non-Lorentzian (table 2). The squared chi for this run is 1.91. When using the Voigt profile as a model and leaving both I-, and r, as free parameters in the fitting procedure, one gets a larger value for r, than the natural width, and for ro a value of 0.229 mm/s (table 2). The squared chi for this run is 1.89. In a second run, where the parameter r, is fixed to its natural value, the extracted r, = 0.292 mm/s is obtained. The value of the squared chi is 1.85, a little smaller than the preceding one. When the broadening spectrum, crossed line in fig. 2c, is fitted using a Gaussian as a model, a value of ro is obtained. The results are shown in fig. 2d and table 2, where one can see that the fitted Gaussian (continuous line) does not completely match the output spectrum of eq. (9). Notice that there is a small puzzling asymmetry in the recovered line. This asymmetry could possibly come from the smoothing operation or it may be an intrinsic feature of this function (see examples in ref. [S]). However, the presence of the pods on the wings makes one think of a Gaussian as the broadening function with a r, as a reported in table 2. According to table 2, the resulting r, resembles the value obtained when using the Voigt profile keeping r, fixed.

4. Conclusions In brief, for values of NT and m larger than those chosen here to give the best results, the FWHMs get

method for small computers

537

narrower than the original ones, and the maximum heights begin to increase in such a way that unrealistic values are obtained. Strong distorted oscillations appear in the experimental recovered spectra although the positions remain unchanged. On the other hand, it is well known that the smoothing procedure causes distortions in the output spectra as other methods do [8,10,17,18]. It is seen, however, that the results obtained with the present method are comparable to those reported in the literature [15], and that as long as the experimental data are not very noisy, the smoothing procedure will not introduce significant distortions in the output spectra; such smoothing procedure may be considered as equivalent to the filtering of the original data [3]. Finally, the removal of the background has shown to be a helpful way for interpreting more satisfactorily the output spectra of eq. (9), where the resulting broadening distribution of the backgroundless experimental data have the same pattern as in the simulated data. A suggested interpretation arising from the above is that the broadening distribution of the Mossbauer spectrum of stainless steel may have a Gaussian shape, where its full width might be in the vicinity of 0.293 mm/s.

References PI M. Migliexini, Nucl. lnstr. and Meth. B36 (1989) 475. 121 Ch. P. Poole Jr., Electron Spin Resonance, 2nd. ed. (Wi-

ley, 1983) p. 459. [31 M.C.D. Ure and P.A. Flinn, Mossbauer

Methodology, Vol. 7 (Plenum, New York, 1971) p. 245. [41 J. Law and D. Hogan, Nucl. Instr. and Meth. B5 (1984) 67. 151 E. SjDntoft, Nucl. Instr. and Meth. 163 (1979) 519. 161 R. Verma, Nucl. Instr. and Meth. 212 (1983) 323. 171 R.A. Brand and G. Le Ca&r, Nucl. Instr. and Meth. B34 (1988) 272. PI I. Vincze, Nucl. Instr. and Meth. 199 (1982) 247. [91 H. Keller, J. Appl. Phys. 52 (1981) 5268. [lOI J. Hesse and A. Rtibartsch, J. Phys. E7 (1974) 526. PII E.O. Brigham, The Fast Fourier Transform (Prentice-Hall, 1974). WI C.J. Batty, SD. Hoath and B.L. Roberts, Nucl. Instr. and Meth. 137 (1976) 179. [13] H. Flores-Llamas, A Cabral-Prieto, H. JimenezDominguez and A. Bravo-Ortega, Nucl. Instr. and Meth. A287 (1990) 557. H. Flores-Llamas, A. Cabral[I41 H. Jimtnez-Dominguez, Prieto and A. Bravo-Ortega, Nucl. Instr. and Meth. A279 (1989) 625. P51 W. Meisel, Phys. Status Solidi B43 (1971) K129. WI Charles S. Duris, ACM Trans. Math. Software 6 (1) (1980) 92. [I71 A.F. Rex, M.A. Olson and N.J. Silva, Appl. Spectrosc. 42 (5) (1988) 775. WV S. Nikolov and K. Kantchev, Nucl. Instr. and Meth. A256 (1987) 160.