Transcript Sophisticated high-resolution deconvolution algorithms for

Sophisticated high-resolution
deconvolution algorithms for
spectroscopic data
V. Matoušek, M. Morháč
Institute of Physics, Slovak Academy of
Sciences, Slovakia
ACAT 2013, Beijing
May 16 – 21, 2013
Introduction
• One of the most delicate problems of any spectrometric method are
those related to the extraction of the correct information out of the
experimental spectra.
• The quality of the analysis of spectrometric data consists in correct
identification of the existence of peaks and subsequently in good
estimation of their positions and intensities (areas).
• The peaks as the main carrier of spectrometric information are very
frequently positioned close to each other.
• The extraction of the correct information out of the spectra sections,
where due to the limited resolution of the equipment, signals coming
from various sources are overlapping, is a very complicated
problem.
• The deconvolution methods represent an efficient tool to improve
the resolution in the data. It is of great importance mainly in the
tasks connected with decomposition of overlapped peaks
(multiplets) and subsequently for the determination of positions and
intensities of peaks in gamma-ray spectra.
• Stationary discrete system that satisfies the superposition principle can be
described by convolution sum
i
y  i    x  k  h  i  k   n(i )  x  i   h  i   n(i ) i  0,1,..., N  1
k 0
where x(i) is the input into the system, h(i) is its impulse function (response),
y(i) is the output from the system, n(i) is additive noise and “ * “ denotes the
operation of the convolution.
In matrix form it can be written
y  Hx  n
where the marix H has dimensions N × M, the vectors y, n have length N
and vector x has length M, while N ≥ M (overdetermined system).
• The problem of finding x, when H, y are known is a discrete ill-posed
problem and requires regularization techniques to get adequate solution.
Three types of regularization methods are very often used:
• smoothing,
• constraints imposition (for example only non-negative data are accepted),
• choice of a prior information probability function - Bayesian statistical approach.
Tikhonov-Miller regularization. The functional
Hx  y
2
 Qx
2
Q, α are the regularization matrix and parameter, respectively,
is minimized. The solution can be obtained by solving the equation

x  H H  Q Q
T
T

1
HT y
Zero-th order or Tikhonov regularization

x  HT H  

1
HT y
Riley Algorithm (commonly called iterated Tikhonov regularization). To obtain smoother
solution one may use algorithm of Tikhonov-Miller regularization with iterative refinement
x
( n 1)
x
(n)

 H H  Q Q
T
T
 H
1
T
y  HT Hx ( n )

x (0)  0
However, because of its iterative nature, the Riley algorithm lends itself to another type of
regularization, so called Projections On Convex Sets – POCS. There exist many regularization
methods based on this kind of regularization. It means we set all negative elements to zero after
each iteration.
Fig.: Solution
obtained
using
Tikhonov method of regularization.
Fig.: Illustration of Riley algorithm of
deconvolution - thick line is original
spectrum, thin line represents
spectrum after deconvolution.
Fig. 64 Riley
deconvolution
with
POCS regularization - thick line is
original spectrum, thin line represents
spectrum after deconvolution.
Van Cittert Algorithm. The basic form of Van Cittert algorithm for discrete convolution
system is



x ( n 1)  x ( n )   H T HH T y  H T HH T Hx ( n )  x ( n )   y'  Ax ( n )
where:
A is system Toeplitz matrix,
n represents iteration number and
μ is the relaxation factor.

Fig. 65 Original
spectrum
(thick
line), deconvolved spectrum using Van
Cittert
algorithm
(without
regularization)
and
deconvolved
spectrum using Van Cittert algorithm
and regularized via POCS method.
Gold Algorithm. If we choose the local variable relaxation factor
and we substitute it into Van Cittert formula we get the
Gold deconvolution algorithm:
x
( n 1)
 i   M 1
y ' i 
A
m 0
im
x
(n )
m
x(n)  i 

i 
x(n)  i 
M 1
A
m 0
im
x(n)  m 

x (0)  1,1, ,1
T
Its solution is always positive when the input data are positive, which makes the
algorithm suitable for the use for naturally positive definite data, i.e., spectroscopic
data.
Fig.: Original spectrum (thick line) and
deconvolved spectrum using Gold
algorithm (thin line) after 10000
iterations. Channels are shown as small
circles.
Fig.: Original spectrum (thick line) and
deconvolved spectrum using Gold
algorithm
(thin line)
after 50000
iterations.
Boosted deconvolution algorithm
• Iterative positive definite deconvolutions (Gold, Richardson-Lucy, Muller and
MAP) converge to stable states.
• It is useless to increase the number of iterations, the result obtained practically
does not change.
• Instead of it we can stop iterations, apply a boosting operation and repeat this
procedure. Boosting operation should decrease sigma of peaks.
• Then, the algorithm of boosted Gold or other iterative deconvolutions is as
follows:
x    1,1,...,1
0
1. Set the initial solution
T
2. Set required number of repetitions R and iterations L
3. Set the number of repetitions r=1
4. According to either Gold, or other iterative algorithms for
n  0,1,..., L  1 find solution x(L)
5. If r=R stop the calculation, else
 0

i   x
a. apply boosting operation, i.e., set x
p is boosting coefficient >0, p  1.1  1.2
b. r=r+1
c. continue in 4.
 L
 i 
p
i  0,1,..., N  1
Fig.: Original spectrum (thick line) and
deconvolved spectrum using boosted Gold
algorithm (thin line) after 200 iterations and 50
repetitions (p=1.2).
Fig.: Original spectrum (thick line) and
deconvolved spectrum using boosted
Richardson-Lucy algorithm (thin line)
after 200 iterations and 50 repetitions.
Fig.: Original spectrum (thick line) and
deconvolved spectrum using boosted MAP
algorithm (thin line) after 200 iterations and
50 repetitions.
Fig. Original gamma-ray spectrum (thick line) and deconvolved spectrum using
classic Gold algorithm (thin line) after 10 000 iterations.
Fig.: Gamma-ray spectrum (thick line) and deconvolved spectrum using boosted
Gold algorithm (thin line) after 200 iterations and 50 repetitions (p=1.2)
Fig.: Experimental gamma-gamma-ray
spectrum (after background elimination).
Fig.: Spectrum after boosted Gold
deconvolution (50 iterations repeated 20
times).
Blind Deconvolution. Sometimes we don’t know the response function, we only
suppose its shape or, respectively, we know theoretically its analytical form, which can
substantially differ from the reality.
We can observe that the operation of the convolution is commutative, which gives us the
possibility of implementation of iterative algorithm of blind deconvolution scheme:
1. Estimate an initial response function h(0)(i) from y(i). This can be accomplished in
different ways. For example, we can choose the narrowest line lobe of y(i), and by
Gaussian fitting of the selected part we can get an initial estimate of the response
function.
2. Accomplish the deconvolution using the estimated response function h(0)(i) and given
y(i)
x(n) i   y i  h(n) i 
where ** denotes the operation of deconvolution.
3. Accomplish the deconvolution using the estimated spectral function for finding an
improved estimation of the response function
h(n1) i   y i  x(n) i 
4. Repeat steps 2, 3 for improving the estimation of x(i), h(i). One can use any of the
presented deconvolution algorithms.
Fig.: Original spectrum (thick line)
and deconvolved spectrum using
blind Gold algorithm (thin line) after
1 000 iterations and 50 repetitions.
Fig.: Deconvolved response function.
Linear systems with changing response
For time dependent linear system analogously as in the case of invariant systems it holds
i
y  i    x  k  h  i , k   n(i )  x  i   h  i   n(i )
i  0,1,..., N  1
k 0
In this case, the columns of the system matrix are not simply shifted as it was in the previous
examples.
In these systems, the shape of response changes with independent variable (time, in the case
of nuclear spectra energy, etc.).
Problem of identification of the system matrix.
Fig.: Principle of interpolation of
response function (EGSnrc simulation
package).
Fig.: Illustration of interpolation between two
simulated responses
Fig.: Synthetized response matrix for
continuum
decomposition
for
experimental data from
Gammasphere.
Fig.: Detail of the response matrix.
Fig.: Measured spectrum of Co56.
Fig.: Spectrum of Co56 after Gold
deconvolution using above displayed
response matrix.
Fig.: Detail of spectrum before deconvolution
(dotted line), after classic (thin line) and boosted
Gold deconvolution (bars).
Conclusions
• We have discussed and analyzed a series of deconvolution methods. There
exists immense number of variations of deconvolution algorithms applied in
various scientific fields, e.g. astronomy, image processing, geophysics, optics
etc. We have introduced only basic classes of deconvolution algorithms that
can be easily employed and implemented for the processing of spectrometric
data.
• We have restricted our investigations only to positive definite methods.
Though they improve substantially the resolution in the spectra they are not
efficient enough to decompose closely positioned peaks.
• We have presented boosted deconvolution algorithms, wich are able to
decompose the overlapped peaks practically to δ functions while
concentrating the peak areas to one channel.
• Though the procedure is fully automatic, due to large variability of the data,
some intervention of the user and tuning of some parameters are required.
• Several algorithms were also implemented in ROOT system in the form of
TSpectrum, TSpectrum2 and TSpectrum3 classes, developed in collaboration
with CERN.
Some relevant publications:
1. M. Morháč , V. Matoušek: Multidimensional Experimental Data Processing
in Nuclear Physics, In:Computer Physics, Eds.: S. Doherty and A. Molloy,
Book Series, Nova Science Publishers, Inc., ISBN: 978-1-61324-790-7,
2012, pp. 1 - 236.
2. M. Morháč, V. Matoušek, High-resolution boosted deconvolution of
spectroscopic data, Journal of Computational and Applied Mathematics,
235 (2011) 1629–1640.
3. M. Morháč, V. Matoušek: Complete positive deconvolution of
spectrometric data. Digital Signal Processing, Volume 19, Issue 3, May
2009, Pages 372-392.
4. M. Morháč, V. Matoušek: Applied Spectroscopy 62, 91 (2008).
5. M. Morháč: Deconvolution methods and their applications in the analysis
of gamma-ray spectra, Nucl. Instrum. Methods Phys. Res., Sect. A 559 (1)
(2006), pp. 119-123.