Application of digital spectral analysis to decomposition of positron lifetime spectra

New method of positron annihilation lifetime spectra processing is proposed. It is based on digital spectral analysis, namely, on Prony method. The regularization procedure to solve the inverse problem of finding parameters of useful signals present in the data is formulated. A new approach to det...

Повний опис

Збережено в:

Бібліографічні деталі
Дата:	2005
Автор:	Kul'ment'ev, A.I.
Формат:	Стаття
Мова:	English
Опубліковано:	Національний науковий центр «Харківський фізико-технічний інститут» НАН України 2005
Назва видання:	Вопросы атомной науки и техники
Теми:	Экспериментальные методы и обработка даных
Онлайн доступ:	http://dspace.nbuv.gov.ua/handle/123456789/81237
Теги:	Додати тег Немає тегів, Будьте першим, хто поставить тег для цього запису!
Назва журналу:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Цитувати:	Application of digital spectral analysis to decomposition of positron lifetime spectra / A.I. Kul'ment'ev // Вопросы атомной науки и техники. — 2005. — № 6. — С. 93-100. — Бібліогр.: 11 назв. — англ.

Репозитарії

Digital Library of Periodicals of National Academy of Sciences of Ukraine

id	irk-123456789-81237
record_format	dspace
spelling	irk-123456789-812372015-05-14T03:02:37Z Application of digital spectral analysis to decomposition of positron lifetime spectra Kul'ment'ev, A.I. Экспериментальные методы и обработка даных New method of positron annihilation lifetime spectra processing is proposed. It is based on digital spectral analysis, namely, on Prony method. The regularization procedure to solve the inverse problem of finding parameters of useful signals present in the data is formulated. A new approach to determine the number of components in a spectrum is proposed. Results on simulated spectra demonstrate the high-resolution capability of the method. Запропоновано новий метод чисельного аналізу спектрів часу життя позитронів. Він базується на цифровому спектральному аналізі, а саме на методі Проні. Сформульовано регуляризаційну процедуру вирішення зворотної задачі знаходження параметрів корисного сигналу, який присутній у даних. Запропоновано новий підхід до визначення числа компонент у спектрі. Результати, отриманні на модельних спектрах, демонструють високу роздільну здатність метода. Предложен новый метод численного анализа спектров времени жизни позитронов. Он основан на цифровом спектральном анализе, конкретно, на методе Прони. Сформулирована регуляризационная процедура решения обратной задачи нахождения параметров присутствующего в данных полезного сигнала. Предложен новый подход определения числа компонент в спектре. Результаты, полученные на модельных спектрах, демонстрируют высокую разрешающую способность метода. 2005 Article Application of digital spectral analysis to decomposition of positron lifetime spectra / A.I. Kul'ment'ev // Вопросы атомной науки и техники. — 2005. — № 6. — С. 93-100. — Бібліогр.: 11 назв. — англ. 1562-6016 PACS: 78.70.Bj http://dspace.nbuv.gov.ua/handle/123456789/81237 en Вопросы атомной науки и техники Національний науковий центр «Харківський фізико-технічний інститут» НАН України
institution	Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection	DSpace DC
language	English
topic	Экспериментальные методы и обработка даных Экспериментальные методы и обработка даных
spellingShingle	Экспериментальные методы и обработка даных Экспериментальные методы и обработка даных Kul'ment'ev, A.I. Application of digital spectral analysis to decomposition of positron lifetime spectra Вопросы атомной науки и техники
description	New method of positron annihilation lifetime spectra processing is proposed. It is based on digital spectral analysis, namely, on Prony method. The regularization procedure to solve the inverse problem of finding parameters of useful signals present in the data is formulated. A new approach to determine the number of components in a spectrum is proposed. Results on simulated spectra demonstrate the high-resolution capability of the method.
format	Article
author	Kul'ment'ev, A.I.
author_facet	Kul'ment'ev, A.I.
author_sort	Kul'ment'ev, A.I.
title	Application of digital spectral analysis to decomposition of positron lifetime spectra
title_short	Application of digital spectral analysis to decomposition of positron lifetime spectra
title_full	Application of digital spectral analysis to decomposition of positron lifetime spectra
title_fullStr	Application of digital spectral analysis to decomposition of positron lifetime spectra
title_full_unstemmed	Application of digital spectral analysis to decomposition of positron lifetime spectra
title_sort	application of digital spectral analysis to decomposition of positron lifetime spectra
publisher	Національний науковий центр «Харківський фізико-технічний інститут» НАН України
publishDate	2005
topic_facet	Экспериментальные методы и обработка даных
url	http://dspace.nbuv.gov.ua/handle/123456789/81237
citation_txt	Application of digital spectral analysis to decomposition of positron lifetime spectra / A.I. Kul'ment'ev // Вопросы атомной науки и техники. — 2005. — № 6. — С. 93-100. — Бібліогр.: 11 назв. — англ.
series	Вопросы атомной науки и техники
work_keys_str_mv	AT kulmentevai applicationofdigitalspectralanalysistodecompositionofpositronlifetimespectra
first_indexed	2025-07-06T05:42:58Z
last_indexed	2025-07-06T05:42:58Z
_version_	1836875076167794688
fulltext	APPLICATION OF DIGITAL SPECTRAL ANALYSIS TO DECOMPOSITION OF POSITRON LIFETIME SPECTRA A.I. Kul'ment'ev Institute of Applied Physics, National Academy of Sciences of the Ukraine, Sumy, Ukraine e-mail: kulmentev@ipfcentr.sumy.ua New method of positron annihilation lifetime spectra processing is proposed. It is based on digital spectral analy- sis, namely, on Prony method. The regularization procedure to solve the inverse problem of finding parameters of useful signals present in the data is formulated. A new approach to determine the number of components in a spec- trum is proposed. Results on simulated spectra demonstrate the high-resolution capability of the method. PACS: 78.70.Bj 1. INTRODUCTION Positron lifetime spectroscopy is one of the best methods for detecting open-volume defects of atomic di- mensions. Positrons implanted into the solid rapidly lose energy via a variety of inelastic processes. The thermal- ized positrons diffuse in the lattice and finally will be annihilated by an electron. The rest energy of electron- positron pair is converted into energy as two annihila- tion γ-quanta with energy 0.511 MeV. Detection of any of these γ-quanta is the signal of positron death. The sig- nal of positron birth depends on the source used but in any case the time difference between these two signals determines the positron lifetime for a single event. The lifetime is a random quantity and thus must be described by its probability distribution function. Many open-volume defects in metals and semicon- ductors are negatively charged and hence form a local attractive potential for positron. Therefore, in the pro- cess of a random walk a positron can be trapped at these defects in a localized bound state. The annihilation rate is determined by the overlapping of positron and elec- tron densities. Negatively charged or neutral open-vol- ume defects have at the same time locally reduced elec- tron density. That is why the lifetime of the positron trapped at these defects increases. In other words, we expect that for open-volume de- fects of k different types the probability distribution function s(t) of the positron lifetime t is a weighted sum of (k + 1) components. It was shown [1] that under quite general assumptions this function has a form 1 1 1( ) exp k j j jj ts t I τ τ + = ж ц = −Ч з чз чи ш е . (1) The aim of the spectrum analysis is the extraction from the experimental signal the set of parameters {τj, Ij} for j = 1,…, k + 1. This analysis allows us to separate out various atomic defect configurations and their relative abundance. The experimental spectrum differs from the function s(t) in several aspects. 1) Function s(t) is the input signal for detecting sys- tem. To obtain the output signal we have to con- volve s(t) with impulse response characteristic of the system. 2) Output signal is formed with the help of a multi- channel analyzer (MCA), that is, the continuous nature of the signal s(t) is discretized. This trans- formation has several consequences. It is a common experience that only one to four com- ponents can be resolved, whereas the number of chan- nels is of order of 100. It means that the problem of spectrum analysis is overdetermined and its solution is usually defined in a least square sense. Events detected in any particular channel form the flow of random events, in the simplest case the Poisson flow. Therefore the spectrum as a whole can be subdi- vided into useful signal and Poisson noise. Discretization of a continuous signal by a multichan- nel analyzer is equivalent (see below) to the procedure of taking sampling values of this signal at some grid points. This grid depends on the signal under study and its properties should be taken into account in spectrum processing. Various approaches have been used to analyze positron lifetime spectra but up to now the most com- monly used are approaches based on the Gauss-Newton non-linear least square fitting [2]. We propose a new ap- proach to the problem and within this approach we stud- ied two problems. The first concerns the determination of the number of components. The conventional solution is to start with one-component fit and add components as long as the variance of the fit decreases. The weakness of this solution is well known in spectral analysis [3]: the vari- ance of the fit usually does not reach a minimum value PROBLEMS OF ATOMIC SCIENCE AND TECHNOLOGY. 2005, № 6. Series: Nuclear Physics Investigations (45), p. 93-100. 93 but rather monotonically decreases as the order of the model is increased. Another approach based on optimal linear filtering and the method of maximum entropy was proposed in [4] but its resolution capability is not high enough and this leads to certain difficulties in treating a spectrum with relatively close components. Another traditional problem concerns the ill-posed nature of the inverse problem [5] under study. This means that the problem is inherently unstable. 2. SAMPLING Let the spectrum reside on a multichannel analyzer with a time increment ∆ per channel. In the experiment the number of counts versus channel number is mea- sured. Therefore we can calculate the total numbers of counts and then for any particular channel [a, b] the probability P of the positron lifetime to fall inside this channel. For the simplest case of a one-component spec- trum this probability is equal to the exponential proba- bility distribution function ( ) ts t Aeα= (2) integrated over t from a to b, ( )( ) . b b a a AP s t dt e eα α α = = −т (3) According to the mean value theorem ( ) ( )( ) ( ), b c a P s t dt s c b a Ae b aα= = − = −т (4) where c ∈ [a, b]. From equations (3), (4) we obtain 1 1ln ec a α α α ∆ж ц−= + з чз ч∆и ш . (5) The distance 1 1ln ec a α δ α α ∆ж ц−= − = з чз ч∆и ш (6) is the same for all channels. It means that the procedure which measures the probability P for the probability dis- tribution function s(t) (2) gives sampled values of this function on an equidistant grid with the same step ∆. Note that the exact position of the grid on the time axis depends on the coefficient α which is usually not known in advance. Therefore, for the positron lifetime spectrum of the form (1) this procedure will determine the sampled val- ues for different components on different grids. Let us suppose that in the model sampled values for the signal s(t) (1) as a whole are defined on the equidistant grid with step ∆ which for simplicity is taken as 0, ,..., ( 1)N∆ − ∆ , (7) where N is the number of channels in the spectrum. It can be easily shown that for a one component sig- nal s(t) (2) two estimators 1 2ˆ ˆ,α α of the damping coeffi- cient α will be the same for two grids: (7) and 0 + δ, ∆ + δ, …, (N – 1) ∆ + δ accordingly. But the estimator 2Â for the amplitude A we have to increase by a factor e-αδ (α < 0) in order to consider formally the estimators 1 2ˆ ˆα α= , (8) 2 2 ˆ 2 1 2 2 ˆ ˆˆ ˆ ˆ 1 A A e A e α δ α α− ∆ ∆ж ц= = з ч−и ш , (9) as obtained for the sampled values of the function s(t) (2) taken at the grid points (7). Similar transformations can be applied to all components in the positron lifetime spectrum (1). 3. COMPUTATIONAL PROCEDURE Experimental positron lifetime spectrum in any chan- nel contains not only the useful signal but also a noise. Therefore let us suppose that the spectrum is defined by the set of sampled values x[1], x[2], …, x[N], taken on the grid (7) and try to estimate this spectrum with a gen- eral p-member model of complex exponentials [ ] 1 ˆ exp ( 2 )( 1) . p j j j j j x n A i f n iα π θ = й щ= + − ∆ +л ые (10) Hereafter, in keeping with the tradition of the signal pro- cessing literature, we use brackets around a discrete variable. In Eq. (10) i is the imaginary unit, 1 ≤ n ≤ N, Aj, αj, fj and θj are amplitude, damping coefficient, frequency and phase constant accordingly for j-th complex expo- nential. Note that we do not impose any restrictions on the order of model p and values of all these parameters. As a consequence sum (10) may present both pure de- creasing exponentials with zero frequencies and un- damped sines with zero damping coefficients. Terms of the first type are the same as the components of the use- ful signal (1) whereas the latter describe the noise part of the spectrum. In shorthand notation we can write down the func- tion of discrete time (10) as [ ] 1 1 ˆ , p n j j j x n h z − = = е (11) where exp( ),j j jh A iθ= (12) exp ( 2 ) .j j jz i fα πй щ= + ∆л ы (13) For N sampled values we have to minimize the sum [ ] 2 1 , N n nρ ε = = е (14) where [ ] [ ] [ ] [ ] 1 1 ˆ , p n j j j n x n x n x n h zε − = = − = − е (15) simultaneously over parameters hj, zj and the number of terms p. It is well known that this is an extremely com- plicated nonlinear problem even for the case where the order p of the exponential model is known in advance. The traditional way to resolve this problem implies the application of iterative algorithms of multidimen- sional minimization. This approach has the following shortcomings. These algorithms are tedious, the ob- tained solution is very sensitive to the choice of starting values of the independent variables and the solution may converge to a local but not to the global extremum in 94 multidimensional space. However these methods as a rule are very flexible and can be applied to any non-lin- ear approximation of sampled values and not necessarily to the exponential one. On the other hand for the latter approximation there is a group of specially developed methods usually defined as the Prony method [6]. The Prony method is not iterative and reduces the non-linear aspects of the problem to finding roots of polynomials. The key point of the Prony method is to consider ex- pression (11) as a general solution of some (yet un- known) linear difference equation [ ] [ ] 0 0, p m a m x n m = − =е (16) with constant coefficients a[m] such that a[0] = 1. Ac- cording to the general theory of these equations the complex constants zj in Eq. (11) are the roots of charac- teristic equation [ ] 0 0. p p m m a m z − = =е (17) Thus, the first stage of computational procedure is to de- termine the coefficients a[m]. There are two distinct cas- es here. If the number of sampled values x[1], …, x[2p] is equal to the number of unknown exponential parame- ters h1, …, hp, z1, …, zp we have a well determined prob- lem and are able exactly to determine them from the set of linear equations (16) writing down for p +1 ≤ n ≤ 2p. In the overdetermined case where N > 2 Eq. (16) should be rewritten in the form [ ] [ ] [ ] 1 , p m a m x n m e n = − =е (18) where e[n] is the error of the forward prediction if the coefficients a[m] are treated as forward prediction coef- ficients. In this case we have to minimize the sum of the squared errors e[n] over all coefficients a[m] so that the latter can be derived as a solution of covariance equa- tions of forward prediction. In the second stage obtained coefficients a[m] are used to form a polynomial (17). Roots of this polynomi- al give the following estimators of damping coefficient αj and frequency fj for j the exponential ln j j z α = ∆ , Im1 . 2 Re j j j z f arctg zπ = ∆ (19) If the complex constants zj are known then Eq. (11) becomes linear with respect to the complex exponential parameters hj. Minimization of sum of the squared errors of the exponential approximation (11) over these param- eters leads to a standard normal equation of least square method ( ) ,H H=Z Z h Z x (20) where superscript "H" denotes Hermitian conjugation, [ ] [ ] [ ] 1 1 2 2 1 1 1 1 2 1 1 1 1 2 , , p N N N pp h x z z z h x h x Nz z z− − − ж ц ж цж ц з ч з чз ч з ч з чз ч= = =з ч з чз ч з ч з чз чз ч з чз ч и ш и ши ш Z h x K K M M O M M M K . (21) The solution of the set of linear equations (20) is the third stage of the computation procedure. Obtained com- plex parameters are then used for estimations of ampli- tude Aj and phase constant θj of j the exponential (see Eq. (12)) j jA h= , Im Re j j j h arctg h θ = . (22) 4. ALGORITHMS At the first stage of computational procedure we have to solve the covariance equations for linear predic- tion. For this purpose we used the fast algorithm initially proposed in [7] and modified later in [8]. This algorithm resembles the Levinson algorithm for Toeplitz matrices but is much more complicated because a Hermitian ma- trix of covariance normal equations is not a pure Toeplitz matrix but can only be represented as a product of matrices of this type. This algorithm like the Levin- son one has recursive structure and therefore for a given order p allows us to obtain the whole set of least square solutions for all orders less than or equal to p. More- over, the algorithm gives not only the solution of covari- ance equations for forward prediction but at the same time the similar solution for backward prediction. At the second stage we used the Jenkins-Traub method [9] for finding the roots of a polynomial. This method finds roots (one per time), roughly speaking in ascending order of their modulus, then deflates found complex root and decreases the order of polynomial. At the third stage we have to solve the set of normal equations of the linear least square method. Even for real coefficients of linear prediction the roots of the characteristic Eq. (17) in the general case will be com- plex and hence we have to solve the set of normal equa- tions with a Hermitian matrix. We used for this purpose the standard algorithm of Cholesky decomposition [10]. 5. RESULTS To test the performance of the program Q_fit we used several series of simulated spectra. This approach is similar to the one proposed in [4]. Moreover, in order to compare the results of processing of the spectrum with different programs we used the same set of input parameters as in [4]. For a given set of lifetime values and relative intensi- ties a number of different spectra were simulated which formed the series. Each spectrum in the series consists of a given number of counts. For each count we first se- lected a state in (1) and then for this state generated one deviate with an exponential distribution function. In both cases the inverse function method was used to transform the random deviate with a uniform probability distribution. Each spectrum has been convoluted with a spectrom- eter resolution function which as in [4] was taken in the form of the single Gaussian with FWHM = 270 ps. For convolution we used a fast Fourier transform algorithm after prior setting up of buffer zone of zero-padded val- ues at the end of a spectrum and bounding the Fourier transform of 95 spectrometer resolution function. Similar algorithms were used for the inverse operation of deconvolution. (The details can be found in [11]). Series A. 2 components: 150 ps 50%, 250 ps 50%; number of channels 128; time calibration 33 ps/channel; counts per spectrum 2×106; number of spectra 5. Let us use this simplest two-component spectrum with well separated and equally intensive components to illustrate the procedure of determining the number of states in (1). First of all we applied the program Q_fit to the pure noise-free spectrum when the sampled values were obtained by exact integration of the useful signal (1) over the channels of MCA. For a successively in- creasing order of exponential model (10) we obtained: 1) p = 1: pure decreasing real exponential (that is both frequency and phase constant are zero within the computational accuracy) and average lifetime τ ≈ 195 ps. 2) p = 2: the single exponential splits into two pure de-creasing real exponentials whose parame- ters are exactly equal to the input ones (τ1 = 150 ps, I1 = 0.5; τ2 = 250 ps, I2 = 0.5). 3) p = 3, 4, 5: all terms of the exponential mod- el can be subdivided into two groups. The first group is comprised of two real exponentials, the same as in 2). The second group includes the re- maining exponentials of the model and all of them within the computational accuracy have zero inten- sity. 4) p = 6: interruption in the subroutine for the fast solution of the covariance equations of linear prediction because in the iterative process the vari- ance of forward prediction became negative. Fig. 1. Application of the program Q_fit to the noise-free spectrum of series A: (a) mean lifetime t of pure decreasing real exponentials with nonvanishing energy vs. order p of model; (b) variance of forward (solid line) and backward (dotted line) prediction In Fig. 1,a mean lifetime values of real exponentials versus the order of exponential model are shown. From this figure it is clear that the spectrum has two compo- nents. The same conclusion can be drawn from the vari- ance of linear prediction versus the order of exponential model dependency (Fig. 1,b). Indeed, for the value p = 2 divides into two regions with sharply different diminish- ing rate of variance of both the forward and backward prediction. That is although for p = 2 the variance, strictly speaking, does not reach a minimum we never- theless are able to figure out when further addition of exponential terms does not lead to any practically suffi- cient improvement of the approximation. Unfortunately, for simulated spectra the situation is quite different. In this case the number of counts on each grid point is subjected to the statistical fluctuation, and the spectrum as a whole consists of the useful signal and a noise. It can be seen from Fig. 2 that for such a spec- trum the variance of both forward and backward predic- tion monotonically decrease as the order of the exponen- tial model increases. Thus, we have to use a different ap- proach to determine the number of states in a positron lifetime spectrum. Fig. 2. Variance of linear prediction for one particular data set from series A Fig. 1,a contains a hint of this approach. As was not- ed, processing of positron lifetime spectrum belongs to ill-posed inverse problem. As a rule for any particular inverse problem we can define two positive functionals FA and FB. The first of them measures the agreement of a model to the data while the second one measures some- thing like the "smoothness" of the desired solution. When FA by itself is minimized the agreement becomes very good, but the solution becomes unstable. In another extreme case, minimizing FB by itself gives smooth solu- tion that has nothing at all to do with the measured data. So, the central idea in inverse theory is to minimize the weighted sum of these two functionals. It can be seen from Fig. 2 that increasing of expo- nential model order leads to the improvement of the agreement between the solution and the underlying spec- trum, that is the solution becomes closer and closer to the nonsmooth initial signal. But the solution sought must look like the useful signal (1). So we have to smooth the derived solution. One way to do this is to maintain in it only a few most energetic pure decreasing exponentials. The energy of a discrete sequence x[i] (i = 1, …, N) from a standpoint of signal processing is defined as 2 1 [ ] . N i E x i = = е (23) Therefore let us subdivide terms of exponential ap- proximation into two groups. The first group will be comprised by pure decreasing real exponentials, that is, by terms j satisfying following conditions 0jα < , jf ε< , ,jθ ε< (24) where ε is the machine-dependent precision parameter which in our calculation was taken as 10-5. Another group will contain the remaining exponentials of the model, that is, terms for which any of the conditions (24) has broken down. In Fig. 3,a for a particular data set from series A the mean lifetime t for all exponentials from the first group versus order of the exponential model is shown. It is seen from this figure that the first group is relatively small. For the order of exponential model from 1 to 63 the size of this group falls in the range 1 (p = 1) to 6 (p = 56). Fig. 3. Application of program Q_fit to the same data set from series A as in Fig. 2: (a) mean lifetime of pure decreasing real terms of the exponential approxi- mation (10) whose energy is ≥ 0.005 Es (large cross) < 0.005 Es (small cross), where Es is the energy of the spectrum; (b) mean lifetime of pure decreasing real ex- ponentials with nonvanishing energy. These exponen- tials form a regularized solution Let us restrict our consideration to the spectra, which do not contain very weak components. In this case we can proceed further with a classification of exponential terms in approximation (10). For any spectrum we can calculate its energy Es and then leave in the first group only those exponentials whose energy exceed some pre- scribed fraction (say, 0.005) of Es (Fig. 3,b). Thus far we have classified only individual terms of the exponential model. But there is also a restriction, which should be imposed on the first group of exponen- tials as a whole. Indeed, these exponentials will form the useful signal (1) and as any probability distribution function this signal should be normalized. Therefore we can impose on each useful component the set of follow- ing conditions (master conditions): 1) it is real and pure decreasing (Eq. (24)); 2) its energy is not vanishingly small as compared with the energy of the spectrum; 3) for given order p these components as a whole form a normalized probability distribution function. Let us define the optimal range of the exponential model as a range of p where the number of its terms sat- isfying master conditions and their mean values are vir- tually constant. It can be seen from Fig. 3,b that for se- ries A the optimal range is defined by the following in- equality 10 ≤ p ≤ 63. Thus, we get the conclusion that a particular spectrum from this series consists of two com- ponents and these two components form the smoothed desired solution. What we have just described can be viewed as a way to determine the number of states in a positron lifetime spectrum and at the same time as a regularization proce- dure applied to the inverse problem under study. It is worthwhile to illustrate the resolution capability of the computational procedure. The energy spectral density of the exponential approximation (11) may be defined by numerous means depending on the adopted assumptions concerning the behavior of this approxima- tion beyond the sampled range. It is generally believed that following two-side function [ ] ( ) 1 * 1 , 0; ˆ , 0, p n j j j p n j j j h z n x n h z n = − = м іп п= н п <п о е е (25) provides the best spectral estimations over other such functions. The z-transformation of the function (25) is given by ( ) ( ) ( ) * 1 * 1 * 2 1 1/ ˆ ( ) 1 1/ / p j j j j j j j j z z z X z h z z z z z z − − − = й щ− к ъ= к ъ− + +к ъл ы е . (26) Hence, according to the general recipe, we can calculate the energy spectral density of the exponential approxi- mation 2 exp( 2 ) ˆ( ) ( ) , z i f E f X z π= ∆ = ∆ Ч (27) where f is the real frequency from the range – fc< f< fc and fс = 1/(2∆) is the Nyquist critical frequency. For the same data set as in Fig. 3 and for the highest possible order of exponential approximation (p = 63) the energy spectral density of this approximation is shown in Fig. 4. For comparison in this figure the ener- gy spectral density of a regularized (smooth) solution is also shown. It can be seen that regularization procedure is equivalent to the extraction of noise from the expo- nential approximation. The spectral density in Fig. 4 illustrates the resolu- tion capability of the computational procedure taken in the ordinary sense, i.e. as a capability to distinguish two sine waves with close frequencies. But we are interested in a different definition of this term as the capability of the computational procedure to distinguish two pure de- ceasing exponentials with close damping coefficients. So let us formally treat frequency f as a complex vari- able f = f + iα and consider the energy spectral density (27) as the section of the right-hand side complex func- tion along the real axis in the complex frequency plane. Fig. 4. Energy spectral density of the exponential approximation of data set from series A (solid line) and of the regularized solution (dotted line) A similar section of this function along imaginary axis is shown in Fig. 5,a for the same order of the expo- nential model (p = 63) as in Fig. 4. For simplicity in- stead of the damping coefficient α we chose in Fig. 5 the mean lifetime τ (τ = –1/α) as an equivalent independent variable. Formally the function in Fig. 5,a represents the result of an analytic continuation of the function in Fig. 4 from the real axis into the complex frequency plane. Fig. 5,a clearly demonstrates the high-resolution ca- pability of the computational procedure to distinguish pure decreasing exponentials. (Note that y-axis in Fig. 5, has a logarithmic scale). From our point of view the high resolution power of the computational procedure may be ascribed to the covariance method of linear prediction and to the special character of the basic set used for the exponential approximation. The covariance method does not use a windowing operation and it is commonly adopted that the absence of this operation leads to the increasing resolution power in any spectral method. It can be seen from the Fig. 3,b that in the range 3 ≤ p ≤ 7 estimations of the average value of two pure decreasing exponentials are biased and the bias mono- tonically decreases as the p is increased. That is why we have to use an exponential model with a relatively high order. But in many spectral methods the higher is the or- der of the model the higher is the probability of splitting of the spectral line [3]. This is also true for the proposed computational procedure as can be seen from Fig. 5,b. In Fig5,b the same function as in Fig. 5,a is shown but calculated for a different order of the exponential model. For p = 60 the first group initially contains four pure de- creasing exponentials (Fig. 3,a). One of them with the shortest lifetime has a very low intensity and hence can be eliminated. The three remaining exponentials satisfy the master conditions but this result takes place for a sin- gle value of p in the range 10 ≤ p ≤ 63 and hence can be taken as accidental. That is why we have to consider the whole range of the exponential model order and make a statistically plausible conclusion, that is, conclusion that takes place for most values of p from this range. It should be noted that the fast algorithm for the solution of covariance equations of linear prediction applied for highest possible order p at Fig. 5. Resolution capability of the Q_fit program for two values of the model order and the same data set from series A as in Fig. 3 the same time gives solutions for all orders less than or equal to p without any additional computation. So con- sideration of the whole range of p does not lead to any significant increasing in the time for spectrum processing. Processing of the series of spectra is naturally subdi- vided into two successive steps. First, for one particular spectrum from the series the above-described procedure is applied. Results are the number of states in the positron lifetime spectrum and the optimal range of the exponential model order. At the second step these data serve as input parameters for the computational proce- dure. The program for a given value of p from the opti- mal range determines the parameters for all terms in the exponential approximation (10), it sorts them into two groups and checks up the first group against the master conditions. If both master conditions are satisfied and the number of states in the regularized solution is equal to the one determined at the first step then this solution is accepted. Otherwise the next value of p from the opti- mal range is selected. In Table 1 are listed the results for series A obtained by MELT [4], Posfit [4] and Q_fit. In this table we give the mean values of lifetimes and intensities as well as the standard deviations. Both MELT and Posfit give reli- able but slightly biased estimates. Q_fit gives very good unbiased results with low errors. From our viewpoint this distinction in estimates may be due to the different procedures of noise treatment adopted in these methods. Table 1. Mean values for lifetimes τi and weights Ii (i = 1, 2) found by MELT, Posfit [4] and Q_fit for simu- lated spectra from the series A. The errors represent the standard deviations from the mean values Analysis method t1(dt1), ps I1 (dI1), % t2(dt2), ps I2 (dI2), % Simulation 150 50 250 50 MELT [4] 148.2 (4.0) 46.7 (3.3) 242.4 (4.0) 53.3 (3.3) Posfi [4] 144.3 (5.6) 45.2 (3.6) 244.8 (3.3) 54.8 (3.6) Q_fit 149.5 (3.1) 50.2 (1.8) 250.4 (1.6) 50.0 (2.0) Table 2. Mean values for lifetimes τi and weights Ii (I = 1,…, 4) found by MELT, Posgauss [4] and Q_fit for simulated spectra from the series D. The errors represent the standard deviations from the mean values Analysis method t1(dt1), ps I1(dI1), % t2(dt2), ps I2(dI2), % t3(dt3), ps I3(dI3), % t4(dt4), ps I4(dI4), % Simulation 100 25 250 25 600 25 1000 25 MELT [4] 103.1 (0.6) 25.27 (0.35) 254.9 (4.3) 23.94 (0.36) 559.3 (14.1) 20.92 (0.80) 947.8 (6.4) 29.87 (0.94) Posgauss [4] 100.6 (0.7) 25.31 (0.40) 254.9 (4.8) 25.34 (0.27) 618.4(13.7) 25.89 (0.57) 1012.1 (8.1) 23.47 (1.02) Q_fit 100.2 (1.7) 24.81 (1.32) 245.9 (14.2) 24.59 (0.78) 592.0 (28.0) 25.34 (0.93) 998.9 (12.9) 25.24 (1.78) Fig. 6. Mean lifetime of the exponential model terms satisfying the master conditions for one particu- lar data set from series D The next two series illustrate the ability of the differ- ent methods to resolve components with close mean life- times. Series B. 2 components: 150 ps 50%, 220 ps 50%; number of channels 128; time calibration 33 ps/channel; counts per spectrum 2×106; number of spectra 5. In this series the two components are closer than in series A. All methods are able to resolve the two components and give estimations with the same qualitative characteristics as in series A. The main difference consists in sufficient increasing of errors for all methods. Series C. 2 components: 150 ps 50%, 190 ps 50%; number of channels 128; time calibration 33 ps/channel; counts per spectrum 2 × 106; number of spectra 5. "Here MELT, unable to resolve the two components at the op- timal entropy weight, showed a wide peak with a mean lifetime of ≈ 169 ps" [4]. Posfit also runs into problems for this series. "The result obtained for the two compo- nent fits are not stable: In the first spectrum the two components are reproduced quite well, in two others the result is approximate, while in the remaining two a spu- rious short or long component appears" [4]. Q_fit are still able to resolve the two components and gives unbi- ased estimators but the errors have further increased up to the ± 10 % level on the mean lifetimes and to the ± 20 % level on the intensities found. Series D. 4 components: 100 ps 25%, 250 ps 25%, 600 ps 25%, 1000 ps 25%; number of channels 128; time calibration 58 ps/channel; counts per spectrum 2× 106; number of spectra 10. This is an artificial and Fig. 7. Resolution capability of the Q_fit program for the same data set from series D as in Fig. 6 complicated data set used mainly to test the performance of the programs. The results found by different methods are summarized in Table 2. All methods resolve four components in spectra but the MELT results show a correlation between two longest components [4]. In this respect both Posgauss ("modified version of Posfit" [4]) and Q_fit performs better. As a whole Q_fit gives esti- mators with the smallest bias but with significantly larg- er errors. The scatter in the results of Q_fit may be due to the high resolution capability of this method. In Fig. 6 for a particular data set from series D the mean lifetime for all exponentials satisfying the master conditions is shown. The tree structure in Fig. 6 is simi- lar to the one for noise-free case but is stretched-out over approximately a ten times longer range of the p-ax- is. The whole range of the exponential model order in Fig. 6 can be subdivided into three intervals: 1) 8 ≤ p ≤ 15; 2) 32 ≤ p ≤ 35; 3) 42 ≤ p ≤ 57, where the parameters of pure decreasing energetic expo- nentials remain virtually constant separated by transition ranges. The first of them corresponds to a three compo- nent model, the second to the four components with bi- ased estimations and the third to the four component model with unbiased estimations. So, the more exactly the noise in the spectrum is described the more reliable is the spectrum decomposition. The optimal range of the exponential model order for the series D is much shorter than for series A and ex- tends from 42 to 57. Fig. 7 illustrates the resolution ca- pability of the computational procedure for one particu- lar order of exponential model from the optimal range. It can be seen from this figure that the useful signal con- sists of four components, which are determined with high-resolution capability. CONCLUSION A new approach to analyze the positron lifetime spectra was proposed. This approach is based on a non- linear least square method but differs from the existing ones in that it uses fast algorithms and a more flexible basic set. This allows to use a model function with a large number of components and thus to describe readily and with reasonable accuracy both useful and noise parts of the spectrum. These components can further be separated and this separation is equivalent to a regularization procedure. It should be stressed that in contrast to the common prac- tice we did not use any low pass filter. Therefore, the proposed method is free from the undesired conse- quences of such filtering. As a result a regularized smooth solution gives unbiased estimations of parame- ters of the useful signal with a high-resolution capability. REFERENCES 1. Positrons in Solids. Ed. P. Hautojärvi. Berlin: “Springer”, 1979, 266 p. 2. P. Kirkegaard et al. Positron annihilation. Singapore: “Wordl Scientific”, 1989, p. 642. 3. S.L. Marple, Jr. Digital Spectral Analysis with Applications. Moscow: “Mir”, 1990, 584 p. (in Russian). 4. A. Shukla et al. Analysis of positron lifetime spectra using quantified maximum entropy and a general linear filter // Nucl. Instr. and Meth. 1993, v. A335, p. 310-317. 5. A. Tikhonov, V. Arsenine. Méthodes de Ré- solution de Problémes Mal Posés Moscow: “Mir”, 1976, 288 p. 6. J.L. Lacoume. Nouvelles méthodes d’anal- yse spectral from J.Max: Traitement du Signal et Application aux Mesures Physiques. Paris: “Mas- son”, 1987, 385 p. 7. M. Morf et al. Efficient solution of covari- ance equations for linear prediction // IEEE Trans. Acoust. Speech Signal Process. 1977, v. ASSP-25, p. 429-433. 8. S.L. Marple, Jr. Efficient least squares FIR system identification // IEEE Trans. Acoust. Speech Signal Process. 1981, v. ASSP-29, p .62- 73. 9. A. Ralston, P. Rabinowitz. A First Course in Numerical Analysis. New York: “Graw-Hill”. 1978. §8.9-8.13. 10. V. Voevodin, Ya. Kuznetsov. Matrix and computation. Moscow: “Nauka”, 1984, 320 p (in Russian). 11. W.H. Press et al. Numerical Recipes in Fortran 90. Cambridge: “Cambridge University Press”, 1992, 1486 p. ПРИМЕНЕНИЕ ЦИФРОВОГО СПЕКТРАЛЬНОГО АНАЛИЗА К ЗАДАЧЕ ДЕКОМПОЗИЦИИ СПЕКТРОВ ВРЕМЕНИ ЖИЗНИ ПОЗИТРОНОВ А.И. Кульментьев Предложен новый метод численного анализа спектров времени жизни позитронов. Он основан на цифро- вом спектральном анализе, конкретно, на методе Прони. Сформулирована регуляризационная процедура ре- шения обратной задачи нахождения параметров присутствующего в данных полезного сигнала. Предложен новый подход определения числа компонент в спектре. Результаты, полученные на модельных спектрах, де- монстрируют высокую разрешающую способность метода. ЗАСТОСУВАННЯ ЦИФРОВОГО СПЕКТРАЛЬНОГО АНАЛІЗУ ДО ЗАДАЧІ ДЕКОМПОЗИЦІЇ СПЕКТРІВ ЧАСУ ЖИТТЯ ПОЗИТРОНІВ О.І. Кульментьєв Запропоновано новий метод чисельного аналізу спектрів часу життя позитронів. Він базується на цифровому спектральному аналізі, а саме на методі Проні. Сформульовано регуляризаційну процедуру вирішення зворотної задачі знаходження параметрів корисного сигналу, який присутній у даних. Запропоновано новий підхід до визначення числа компонент у спектрі. Результати, отриманні на модельних спектрах, демонструють високу роздільну здатність метода. PACS: 78.70.Bj 2. SAMPLING 3. COMPUTATIONAL PROCEDURE 4. ALGORITHMS А.И. Кульментьев О.І. Кульментьєв

Application of digital spectral analysis to decomposition of positron lifetime spectra

Репозитарії

Схожі ресурси