- Email: [email protected]
An Application of Subspace-Based Techniques to Nuclear Spectroscopy
Copyright © IFAC System Identification, Kitakyushu, Fukuoka, Japan, 1997
AN APPLICATION OF SUBSPACE-BASED TECHNIQUES TO NUCLEAR SPECTROSCOPY
Sergio Bittanti, Emilio Gatti, Giancarlo Ripamonti, Sergio M. Savaresi
Politecnico di Milano. Dipanimento di Elettronica e InJormazione. Piazza Leonardo da Vinci.32. 20133 Milano. ITALY.
Abstract. In this paper, a method for the identification of the transfer function of an analog amplifier for nuclear spectroscopy is presented. The proposed technique is based upon a Subspace-based System State-Space Identification method, complemented with a pre-processing technique of the data, which improves the accuracy of the estimation. The identification method is tested both on simulated and real data sets. Key Words. Hybrid filters; subspace-based identification; impulse response; nuclear spectroscopy.
required in order to avoid the so-called "pile-up" phenomenon, i.e. the superposition of two successive impulse responses, which would lead to unbearable errors in the measure of the impulse response peak (Bingefors (1993), Ripamonti et at. (1994» .
1. INTRODUCTION AND PROBLEM STATEMENT
In nuclear spectroscopy, an important problem is that of measuring the energy of pulse signals (which are the typical signals associated, for instance, to the decay of a particle). The measurement of the energy of such impulsive signals is commonly made by means of analog filters. To this purpose, the filter is fed with the impulse, and the height of the maximum peak of the associated output in the time domain is measured. Such a height is indeed linearly proportional to the energy of the impulse.
In order to reduce the cost of such filters, and to give more flexibility to their design, an hybrid filter has been recently developed and proposed in Georgiev and Gast (1993) and Ripamonti et at. (1994). Hybrid filters proved to be useful in various instances, e.g. for LCH experiment at CERN (Bingefors (1993» . In hybrid filters, the pre-filter (F(s» is an analog filter characterised by an infinite impulse response, while the "shaping-filter" ( H(z) ) is constituted by a digital processor (preceded by an AID converter).
The analog filters which are resorted to are usually constituted by the cascade of two sub-filters: the first filter (the "pre-filter") is essentially a low-pass filter having the aim - roughly speaking - of "spreading and flattening" the time domain shape of the impulsive signal, while the second filter (the "shaping filter") takes care of carefully shaping the overall impulse response. This shaping essentially consists in the cut of the log "tail" of the impulse response of the pre-filter: such a feature is strongly
In Ripamonti et at. (1994) the design of the filter H(z) is discussed: the main conclusion is that a "well-shaped" finite time impulse response .can be obtained by means of a FIR-like filter H(z) , and that the so-obtained hybrid filter outperforms the classical analog filters, both in terms of flexibility
1591
and q, one cannot appeal to a general theory; however, a wise choice seems to be that of using d ~ q , d being not much larger than q (i.e.
Hqd
"almost" square).
(iii) Compute the Singular Value Decomposition (SVD) of Hqd : (l)
where Q and V are orthogonal matrices qxq and d xd having dimension respectively, and S is a diagonal matrix with the singular values in non-increasing order. (iv) Truncate the matrices Q,
V,
and
S.
In the
ideal case (i.e. when no noise is present),
Hqd
The most attractive feature of the subspace-based method above described is likely to be its numerical stability and efficiency. As a matter of fact, the most critical and computationally expensive calculus of the subspace-based methods is the SVD decomposition of the extended Hankel matrix. As is well known (see e.g. Golub and van Loan (1989)), the computation of the SVD of an high-dimensional matrix is a numerically stable operation, which can be efficiently performed. This feature is extremely appealing, especially if one compares with the standard parametric identification methods (see e.g. Ljung (1987)), which usually require the minimisation of a multi-variable nonlinear loss function, which calls for iterative optimisation methods, so that one encounters non-trivial convergence and numerical stability problems. Another distinctive feature of the 4SID method is that no initialisation, nor a-priori knowledge of the structure of the system to be identified is required. However such a feature might be regarded as a pitfall, when a quite accurate a-priori estimate of the system to be identified is available: in this case, the subspace-based identification methods are not capable of taking advantage of such an a-priori knowledge, so resulting in a "waste" of information.
is of rank n (the order of the system to be identified), so that the only first n singular values are non-zero. In case of noisy measurement, all singular values will be positive with probability 1. The user must then tackle the non-trivial task of deciding the number of "significant" singular values. Once this decision is performed, the matrices of the SVD decomposition (1) can be partitioned into the "signal" (index "S') and "noise (index "N') parts as follows:
On the other hand, probably the most obnoxious drawback of the 4SID methods is their sensitivity to the intensity of the noise corrupting the signal. To this purpose, in the following subsection, a simple but effective procedure is proposed, having the aim of increasing the Signal-to-Noise-Ratio (SNR) of the original measurement of the impulse response, h(k), in the special case when an a-priori e'stimate of the poles-zeros configuration is available.
where Qs and v~ contain the n principal singular vectors, whose correspondent singular values are collected in the matrix Ss' (v) Compute the matrices:
r = Q-sSsy, ' n = SsY,Vs·
2.2. Pre-processing of the data set.
T
The pre-processing procedure we propose relies on the hypothesis that the signal hO is band-limited, and its band is narrower than that of the corrupting noise. The basic idea of the method is that of filtering the difference between the measured impulse response and an a-priori estimate of the impulse response (which is supposed to be available), and successively re-constructing a new "noisy" impulse response, characterised by a higher SNR than the original one.
These matrices are the estimated extended observability and reachability matrices, respectively. (vi) Compute realisation
A=(r
the
estimate
of
a
state-space
{A, B, C} of the system:
( I:q- I . I",)
)tr
(2 :q,I:. ) ,
B=n
(1", .1:1)'
c=r
(1:1.1:.) ,
where the symbol A(a'h,<".d ) indicates the submatrix of A constituted by its rows of indices {a,a+l, oo. ,b-l,b}, and by its columns of
Pre-Processing Procedure With reference to a transfer function F(s) , and to the associated (sampled) impulse response h(k), k = 1,2' 00' M, consider again a noisy measure
indices {c,c+1,00.,d-1,d}, and the symbol
"t" indicates the operation of pseudoinverse.
h(k) = h( k) + e( k)
(vii) Compute the estimated transfer function F(z) :
of its impulse response, where
e(k) is the corrupting noise.
Fez) = C(zI - Ar B , l
Suppose that an estimate, say P(s) , of the transfer function is available, the correspondent (sampled)
from which all poles and zeros can be evaluated.
1593
AN APPLICATION OF SUBSPACE-BASED TECHNIQUES TO NUCLEAR SPECTROSCOPY
Sergio Bittanti, Emilio Gatti, Giancarlo Ripamonti, Sergio M. Savaresi
Politecnico di Milano. Dipanimento di Elettronica e InJormazione. Piazza Leonardo da Vinci.32. 20133 Milano. ITALY.
Abstract. In this paper, a method for the identification of the transfer function of an analog amplifier for nuclear spectroscopy is presented. The proposed technique is based upon a Subspace-based System State-Space Identification method, complemented with a pre-processing technique of the data, which improves the accuracy of the estimation. The identification method is tested both on simulated and real data sets. Key Words. Hybrid filters; subspace-based identification; impulse response; nuclear spectroscopy.
required in order to avoid the so-called "pile-up" phenomenon, i.e. the superposition of two successive impulse responses, which would lead to unbearable errors in the measure of the impulse response peak (Bingefors (1993), Ripamonti et at. (1994» .
1. INTRODUCTION AND PROBLEM STATEMENT
In nuclear spectroscopy, an important problem is that of measuring the energy of pulse signals (which are the typical signals associated, for instance, to the decay of a particle). The measurement of the energy of such impulsive signals is commonly made by means of analog filters. To this purpose, the filter is fed with the impulse, and the height of the maximum peak of the associated output in the time domain is measured. Such a height is indeed linearly proportional to the energy of the impulse.
In order to reduce the cost of such filters, and to give more flexibility to their design, an hybrid filter has been recently developed and proposed in Georgiev and Gast (1993) and Ripamonti et at. (1994). Hybrid filters proved to be useful in various instances, e.g. for LCH experiment at CERN (Bingefors (1993» . In hybrid filters, the pre-filter (F(s» is an analog filter characterised by an infinite impulse response, while the "shaping-filter" ( H(z) ) is constituted by a digital processor (preceded by an AID converter).
The analog filters which are resorted to are usually constituted by the cascade of two sub-filters: the first filter (the "pre-filter") is essentially a low-pass filter having the aim - roughly speaking - of "spreading and flattening" the time domain shape of the impulsive signal, while the second filter (the "shaping filter") takes care of carefully shaping the overall impulse response. This shaping essentially consists in the cut of the log "tail" of the impulse response of the pre-filter: such a feature is strongly
In Ripamonti et at. (1994) the design of the filter H(z) is discussed: the main conclusion is that a "well-shaped" finite time impulse response .can be obtained by means of a FIR-like filter H(z) , and that the so-obtained hybrid filter outperforms the classical analog filters, both in terms of flexibility
1591
and q, one cannot appeal to a general theory; however, a wise choice seems to be that of using d ~ q , d being not much larger than q (i.e.
Hqd
"almost" square).
(iii) Compute the Singular Value Decomposition (SVD) of Hqd : (l)
where Q and V are orthogonal matrices qxq and d xd having dimension respectively, and S is a diagonal matrix with the singular values in non-increasing order. (iv) Truncate the matrices Q,
V,
and
S.
In the
ideal case (i.e. when no noise is present),
Hqd
The most attractive feature of the subspace-based method above described is likely to be its numerical stability and efficiency. As a matter of fact, the most critical and computationally expensive calculus of the subspace-based methods is the SVD decomposition of the extended Hankel matrix. As is well known (see e.g. Golub and van Loan (1989)), the computation of the SVD of an high-dimensional matrix is a numerically stable operation, which can be efficiently performed. This feature is extremely appealing, especially if one compares with the standard parametric identification methods (see e.g. Ljung (1987)), which usually require the minimisation of a multi-variable nonlinear loss function, which calls for iterative optimisation methods, so that one encounters non-trivial convergence and numerical stability problems. Another distinctive feature of the 4SID method is that no initialisation, nor a-priori knowledge of the structure of the system to be identified is required. However such a feature might be regarded as a pitfall, when a quite accurate a-priori estimate of the system to be identified is available: in this case, the subspace-based identification methods are not capable of taking advantage of such an a-priori knowledge, so resulting in a "waste" of information.
is of rank n (the order of the system to be identified), so that the only first n singular values are non-zero. In case of noisy measurement, all singular values will be positive with probability 1. The user must then tackle the non-trivial task of deciding the number of "significant" singular values. Once this decision is performed, the matrices of the SVD decomposition (1) can be partitioned into the "signal" (index "S') and "noise (index "N') parts as follows:
On the other hand, probably the most obnoxious drawback of the 4SID methods is their sensitivity to the intensity of the noise corrupting the signal. To this purpose, in the following subsection, a simple but effective procedure is proposed, having the aim of increasing the Signal-to-Noise-Ratio (SNR) of the original measurement of the impulse response, h(k), in the special case when an a-priori e'stimate of the poles-zeros configuration is available.
where Qs and v~ contain the n principal singular vectors, whose correspondent singular values are collected in the matrix Ss' (v) Compute the matrices:
r = Q-sSsy, ' n = SsY,Vs·
2.2. Pre-processing of the data set.
T
The pre-processing procedure we propose relies on the hypothesis that the signal hO is band-limited, and its band is narrower than that of the corrupting noise. The basic idea of the method is that of filtering the difference between the measured impulse response and an a-priori estimate of the impulse response (which is supposed to be available), and successively re-constructing a new "noisy" impulse response, characterised by a higher SNR than the original one.
These matrices are the estimated extended observability and reachability matrices, respectively. (vi) Compute realisation
A=(r
the
estimate
of
a
state-space
{A, B, C} of the system:
( I:q- I . I",)
)tr
(2 :q,I:. ) ,
B=n
(1", .1:1)'
c=r
(1:1.1:.) ,
where the symbol A(a'h,<".d ) indicates the submatrix of A constituted by its rows of indices {a,a+l, oo. ,b-l,b}, and by its columns of
Pre-Processing Procedure With reference to a transfer function F(s) , and to the associated (sampled) impulse response h(k), k = 1,2' 00' M, consider again a noisy measure
indices {c,c+1,00.,d-1,d}, and the symbol
"t" indicates the operation of pseudoinverse.
h(k) = h( k) + e( k)
(vii) Compute the estimated transfer function F(z) :
of its impulse response, where
e(k) is the corrupting noise.
Fez) = C(zI - Ar B , l
Suppose that an estimate, say P(s) , of the transfer function is available, the correspondent (sampled)
from which all poles and zeros can be evaluated.
1593