- Email: [email protected]
Constrained digital matched filter method for optimum filter synthesis
ARTICLE IN PRESS
Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351 www.elsevier.com/locate/nima
Constrained digital matched filter method for optimum filter synthesis Wen Xiangyang, Wei Yixiang Department of Engineering Physics, Tsinghua University, 212, Beijing, China Received 29 November 2005; received in revised form 19 December 2005; accepted 19 December 2005 Available online 17 January 2006
Abstract We present a new method to directly calculate the optimum filter in presence of any additive stationary noise, with arbitrary time and domain constraints (flat-top, zero-area, etc.). A more concise re-deduction of digital penalized LMS method (DPLMS) is given. This method is fully developed, and synthesis results of a typical situation are given and compared with the DPLMS method. Optimum filter can be synthesized without a prior knowledge of the noise power spectral density, which makes it suitable to be used in adaptive, selfcalibrating digital spectroscopy. r 2006 Elsevier B.V. All rights reserved. PACS: 29.30.Kv; 02.10.Sp Keywords: Digital filter; Digital spectroscopy; Optimum filter synthesis
1. Introduction Digital pulse processing systems have been broadly used in high-resolution spectroscopy. One of the major causes of inaccuracy in such systems is the electronic noise introduced by the setup, and the signal-to-noise ratio (SNR) can be greatly improved by appropriate digital processing [1,2]. So digital filter plays a crucial role on the system performances, therefore how to synthesize the optimum filter is a key problem in the design of digital spectroscopy systems. When referring to synthesizing ‘‘optimum filter’’, firstly, we should answer the following question: what type does the optimum filter belong to in the sense of SNR? Gatti and Redeka had proved that time-invariant linear filter can achieve the maximum SNR even without the constraint of linearity or time-invariant, in other words, if we find a filter which achieves the maximum SNR in the scope of timeinvariant and linear filters, it is the optimum filter in general [3,4]. It is impossible, at least in principle, to find a better filter even in a more broad scope based on the SNR criterion. Thus in the following we find the optimum filter Corresponding author. Tel.: +86 10 6279 4799.
E-mail address: [email protected] (W. Xiangyang). 0168-9002/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2005.12.199
from time-invariant linear filters. In many situations arising in practice, the optimum filter not only yields the best SNR, but also should satisfies additional constraints imposed from time domain or frequency domain, such as, flat-top needed for minimizing ballistic deficit effect, overall zero-area of the pulse response for baseline rejection [5,6]. There are two types of linear digital filter namely finite impulse response (FIR) filter and infinite impulse response filter. FIR filter uses only current and past input samples to obtain current output value, and this characteristic is very important, since it means that the input pulses occurring N T times ago will have no contribute to the current output, where N is the length of FIR filter and T is the sample interval. Also FIR filter is not sensitive to filter coefficient quantization and easy to implementation. So FIR filters are very suited for digital spectroscopy applications, although they need more computing power, which is not a problem for current DSP processors and FPGA. Extensive researches have been carried out in the past years, and many methods are proposed to synthesize continuous-time weighting functions which cannot be used directly, and digital filters are derived by sampling its analogue counterpart which are sub-optimum [7]. On the
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
other hand, direct synthesis digital filters is more straightforward and easier to deal with. The DPLMS method introduced by Gatti et al. [8] is a general procedure to directly find the optimum digital filters, which implements a relaxed constraints fulfillment technique which allows achieving constraint accomplishments up to the required level of precision, to avoid further and needless degradation in measurement resolution. Many methods are summarized and fully critically reviewed in Ref. [8]. If several constraints exist at the same time, filter synthesized by DPLMS cannot perfectly satisfy all the constraints. In some practical applications, that every constraint is perfectly fulfilled is more preferable, we propose a new method to meet this need, and how constraints are imposed on the filter is explained in detail. 2. Constrained digital matched filter method
347
point N is, Ps ðNÞ ¼ jhT sj2 .
(2)
Since noise is stationary, output noise power is the same at every point, Pv ðNÞ ¼ EjhT vðN : 1Þj2 ¼ hT EfvðN : 1ÞvðN : 1ÞT gh ¼ hT Rv h
ð3Þ
where Rv is the noise autocorrelation matrix, 0 1 rð0Þ . . . rðN 1Þ B C .. .. .. C Rv ¼ B . . . @ A rðN 1Þ rð0Þ where rðnÞ is autocorrelation function, Rv is positive definite. The SNR at point N is jhT sj2 hT R v h
Constrained digital matched filter (CDMF) method is based on the idea that the optimum FIR filter should maximize the SNR function, and at the same time, satisfy all the constrains imposed on by linear equations. Thus the search of the optimum filter is equivalent to solving an equality-constrained optimization problem.
where SNRðhÞ is a function of h.
2.1. Derivation of SNR
A FIR filter’s output is the convolution of input signal with its impulse response, the output is,
The input sequence obtained by sampling the output of the preamplifier can be written as,
ysðnÞ ¼ sðnÞ hðnÞ ¼
xðnÞ ¼ sðnÞ þ vðnÞ where sðnÞ is the input signal sequence with known shape, and vðnÞ is input zero-means noise sequence which is uncorrelated with sðnÞ. The input signal shape can be get as a prior knowledge, or by sampling and averaging [8] or by pole-zero identification method [9]. The CDMF method uses only noise autocorrelation function, which can be calculated from the power spectrum and system functions before the ADC, or as in many situations arising in practice that we have no exact knowledge about the noise, can be estimated with Eq. (1) by sampling a long noise sequence. rðmÞ ¼
1 N jmj
N1m X
xðnÞxðn þ mÞ.
(1)
n¼0
The output of the N-taps FIR with impulse response h ¼ ðhð0Þ; . . . ; hðN 1ÞÞ is, yðnÞ ¼ hðnÞ xðnÞ ¼ hT sðn : n N þ 1Þ þ hT vðn : n N þ 1Þ. From the input signal sðnÞ we choose a N-point segment which contains the maximum power among all the N-point long segments. For convenience we suppose that this segment starts from point one, and define s as s ¼ ðsðNÞ; . . . ; sð1ÞÞ. Among the output sequence, point N should have the maximum SNR. The signal power at
SNRðhÞ ¼
(4)
2.2. Time domain constrains
N 1 X
hðkÞsðn kÞ.
(5)
k¼0
From Eq. (5) time domain constraints on the output signal shape can be easily applied. Here we give an example on how flat-top constrain is implemented. If a M point flat-top is required, we force the successive (M+1) output to be equal, that is to say, ysðNÞ ¼ ¼ ysðN MÞ. Thus we get M equations, ysðN M þ 1Þ ¼ ysðN MÞ ysðNÞ ¼ ysðN MÞ. Substitute ysðnÞ with Eq. (5), and rewrite it in matrix form, 0 1 xðN m þ 1Þ xðN mÞ . . . xðm þ 2Þ xðm þ 1Þ B C .. .. .. B C B C . . . @ A xðNÞ xðN mÞ xð1Þ xðm þ 1Þ 1 0 hð0Þ 0 1 0 C B C B .. C B C B . C ¼ B ... C B C @ A B B hðN 1Þ C A @ 0 namely, At h ¼ bt .
(6)
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
348
So the flat-top constraint is converted to constraint of linear equations about h. 2.3. Frequency domain constrains
with the null space of AðLT Þ1 , according to the relation between null space and range space, H is the orthogonal complement of range of L1 AT [10]. We performs the QR decomposition on L1 AT , L1 AT ¼ QR ¼ ðq1 ; . . . ; qm ; qmþ1 ; . . . ; qN Þ R. |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}
The frequency response of a FIR filter is,
T
HðoÞ ¼ DTFTðhðnÞÞ ¼
N1 X
RðA L Þ
hðkÞe
jko
(7)
k¼0
which is fully decided by hðnÞ, thus any frequency domain constraints on the filter are transformed to constraints on the filter coefficients. Here we give an example on how impose constraint on the filter to remove baseline shift. The meaning of removing baseline is that output will be zero if input stay on a constant level, that is to say, the frequency response at 0 Hz is 0, from Eq. (7) we get, N 1 X
hðkÞ ¼ 0
k¼0
rewrite it in matrix form, Af h ¼ bf
(8)
where Af ¼ ð1; . . . ; 1Þ; bf ¼ 0. We combine the foregoing two constraints as follows: Ah ¼ b 0
(9) At
1
B C A ¼ @ A f A;
0
bt
1
The columns of Q can be divided into two groups: the first M columns of Q form an orthonormal basis of RðAT L1 Þ, since RðAT L1 Þ is the orthogonal complement of subspace H, so the last N M columns of Q form an orthonormal basis of H, H ¼ ðqmþ1 ; . . . ; qN Þ. Thus we get the projection matrix onto H, PH ¼ HH T .
(15) L1 v s
onto H, the residual The project matrix PH projects part is orthogonal with any vector in H, thus L1 v s is divided into two parts: one equal to PH L1 s lying in H, the v other is orthogonal with H, therefore is orthogonal with LTv h. We get, T T 1 ðLTv hÞT ðL1 v sÞ ¼ ðLv hÞ ðPH ðLv sÞÞ.
Eq. (11) can be rewritten as, SNRðhÞ ¼
2 jðLTv hÞT PH ðL1 v sÞj ðLTv hÞT ðLTv hÞ
(16)
which, according to the Cauchy–Schwartz inequality, attains its maximum
B C b ¼ @ bf A ¼ 0
2 SNRmax ¼ kPH L1 v sk
when the filter satisfies
where A has full row rank.
(17) LTv h
¼
cPH L1 v s,
h ¼ cðLTv Þ1 PH L1 v s,
2.4. Solve the CDMF problem From Eqs. (4) and (9), we get an equality-constrained maximization problem, namely, hT ssT h max SNRðnÞ ¼ T h2RN h Rv h
(14)
H
1
or equivalently (18)
where c is any constant. The optimum filter can be scaled in any desirable way. 2.5. Remarks on CDMF method
s:t:
Ah ¼ 0.
(10)
which, according to the Cauchy–Schwartz inequality, attains its maximum when the filter satisfies
If the noise is white, Eq. (13) becomes to h ¼ cs, the filter is a scaled time-reverse replica of the input signal shape, this property resulted in the term matched filter. Eq. (13) provides the optimum filter for color additive noise. The optimum filter with constraints is synthesized by three cascaded steps: whitening, projection and matching: whitening step transforms the noise autocorrelation Rv into a unit matrix, and projection step projects the whitened input signal onto the space determined by the constraint equation, and matching is the last step.
LTv h ¼ cL1 v s
3. New deduction of DPLMS method
Perform Cholesky decomposition on Rv , we get, Rv ¼ Lv LTv . Rewrite Eq. (4) as, SNRðhÞ ¼
2 2 jhT Lv L1 jðLTv hÞT ðL1 v sj v sÞj ¼ T T h Lv LTv h ðLTv hÞ ðLTv hÞ
(11)
(12)
where c is any constant. Without the constraints, we get 1 h ¼ cðLTv Þ1 L1 v s ¼ cRv s.
(13)
We define H as the subspace spanned by all possible LTv h, since h is constrained by Eq. (9), L1 v s may not lie in H, in other words, Eq. (12) may have no solution. H is identical
To begin with, the main idea of this method is recalled: the optimum filter is found by minimizing the weighted sum of the output noise variance and other terms that penalize constraint violations, all penalty terms are in quadratic form. In case of only flat-top and baseline
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
rejection constraints are concerned, the problem can be expressed as, 8 T 8 PV ðNÞ ) 0 > h Rv h ) 0 > > > > > T > < < Y ðNÞ ) 1 kh s 1k2 ) 1 ¼) > > A t h ) bt kAt h bt k2 ) 0 > > > > > : : kA h b k2 ) 0: A f h ) bf f f Assume the weights are ½a1 ; a2 ; a3 ; a4 , the objective function can be written as, FðhÞmin ¼ a1 hT Rv h þ a2 khT s 1k2 þ a3 kAt h bt k2 þ a4 kAf h bf k2 ¼ hT Ah 2bT h þ a2 þ a3 bTt bt þ a4 bTf bf T
a3 ATt At
a4 ATf Af ,
ð19Þ T
where A ¼ a1 Rv þ a2 ss þ þ b ¼ a2 s T þ T T a3 bt At þ a4 bf Af . Since Rv is positive definite and symmetric and other terms are real and symmetric, so A is nonsingular, the solution of Eq. (19) is, DFðhÞ ¼ Ah b ¼ 0 ) h ¼ A1 b.
(20)
349
eliminates the need for any flat-top, as is the case of the present example. A special point should be noticed is that, in digital systems, a flat-top further by one sampling time is needed because of the random jitter between sampling comb and pulse arrival time [11]. In the following, a 15 points flat-top constraint is assumed which is not really needed for the present detector, but just to exhibit the capability of flat-top filters synthesizing. The optimum filters for four cases are synthesized respectively, A: no constraint, B: flat-top constraint, C: baseline rejection constraint, D: flat-top and baseline rejection constraint. The synthesis results are plotted in Fig. 1, it can be seen that output A in Fig. 1(a) is more cusp-like than its counterpart in Fig. 1(b), which is in accord with the matched filter theory. 4.3. Comparisons between DPLMS and CDMF It is easy to verify from Eq. (18) that the synthesis result strictly satisfies all the constraints, so the CDMF method
4
x 104 A B C D
4. Synthesis results and discussions 3 2 Output
In the following, the CDMF method described in Section 2 was applied to a typical situation, the DPLMS method re-deduced in Section 3 was also applied for comparison.
1
4.1. Input conditions
0
We take the small silicon detector as an example, this simulated system has the following conditions: d input current pulse, only series noise and parallel noise with noise corner time tc . The sample frequency is 10 MHz, and the analog system before the ADC is the cascade of three parts: charge-sensitive amplifier, whitening differentiator with time constant followed by a simple R-C anti-aliasing filter, the transfer function is, K ðs þ 1=tc Þðs þ 1=ta Þ
pa2 jtj=ta e 2ta
0
10
20
30
(a)
40 50 time(n)
60
70
80
90
8 A B C D
6
(21) 4
where tc ¼ 0:5 ms; ta ¼ 1:0 ms, tc and ta is the time constant of the differentiator and the anti-aliasing respectively. The overall autocorrelation function is, RðtÞ ¼
−2
(22)
where a2 ¼ 1:52e015 v2 s.
2 Output
F ðsÞ ¼
−1
0 −2 −4
4.2. Synthesis results −6
When designing optimum filter for a given system, constraints should be imposed according to the specific conditions, in cases where the charge collection time of the detector can be regarded to as negligible, which in practice
(b)
0
10
20
30
40
50
60
70
80
90
time(n)
Fig. 1. Comparison of the weight functions with different noise situations: (a) color noise described by Eq. (22), (b) white noise.
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
350
can get perfect flat-top and absolutely baseline rejection at the same time. While from Eq. (20), if no constraints exist, the DPLMS method leads to matched filter and gets the same result as CDMF. If one or more constraints are imposed on, the result is a compromise between the SNR and fulfillment of the constraints. If there is only one constraint, by emphasizing the associated weight coefficient the constraint can be almost perfectly satisfied, at the expense of a degraded SNR, which can be seen from Fig. 2(a) while the flat-top weight value changes from less effective ða3 ¼ 30Þ to highly effective ða3 ¼ 1000Þ, the top of the output is more and more flat. However, if there are more than one constraints, the constraints cannot be perfectly satisfied at the same time, which can be seen from Fig. 2(b), the input sequence staying on a dc voltage of 2, is filtered by filters synthesized with different a4, the output shift is 0.2 when a4 ¼ 1, and is almost completely removed when a4 ¼ 10, but while the ability of removing baseline enhanced, the flatness of the output top is worsened.
1.0015
α3 = 30 α3 = 100 α3 = 1000 CDMF
1.001
1.0005
1
Theoretic SNR for several situations is computed for the filters synthesized by the two methods, with the foregoing input conditions, and is normalized to the SNR given by CDMF with no constraint. Some result is given in Table 1. When both the constraints exist, the two methods give almost the same SNR, but, in practice applications, the imperfect fulfillment of the constraints by DPLMS would worsen the SNR, which is also proved by simulation. 4.4. Constraints’ influence on SNR We define the absolute amplitude-to-rms noise ratio as Z ¼ SNR1=2 , from Eq. (17), the possible maximum Z with given input s is, Z1 ¼ limN!1 kL1 v sk, for convenience, we define the relative amplitude-to-rms noise ratio as Zr ¼ Z=Z1 . When synthesizing the optimum filter, the first decision we have to make is how long the FIR filter is, for this purpose we calculate Zr as a function of FIR length N, the result shows that relatively short FIR length is sufficient to achieve SNR fairly close to the ideal SNR Z1 (Zr 40:99 with FIR length 30). Constraints imposed from the time domain or frequency domain reduce the filter freedom degrees, therefore would results in a degraded SNR. The influence on Zr of various filter length and constraints is plotted in Fig. 3. There are six curves in this figure. In fact, because the three curves of filter length 30; 50; 80, respectively, without baseline rejection are almost identical, which are labeled as ‘Non-BLR’ as a whole, we can conclude that without baseline rejection constraint, filter length has little impact on the SNR when Table 1 The SNR of DPLMS with both constraints is 0.626, the SNR of DPLMS is computed for four cases, ½a1 ; a2 is fixed to ½1; 1
0.9995
2
3
4
5
6
7
8
9
½a3 ; a4
½0; 1
½1; 0
½1; 1
½10; 10
SNR
0.876
0.818
0.626
0.626
n
(a) 5
α4 = 1 α4 = 10 α4 = 1000 CDMF
0.4
4
0.2
3
0
1 Non−BLR 30+BLR 50+BLR 80+BLR
0.9
2
0.8
1 0.7
−1
4.6
−2
0.5
4.2 4
0.4
3.8
−4 −5
0.6
4.4
−3
(b)
ηr
0
0.3 0
20
40
60
80
100
120
Fig. 2. Filters synthesized for the following conditions: flat-top width 6, FIR length 60, weight value is 1 if not explicitly asserted: (a) flat-top synthesized with different value of a3, (b) baseline remove effects with different value of a4.
0.2
0
1
2
3
4
5
6 7 8 9 10 11 12 13 14 15 Flat top Length
Fig. 3. Filter length and constraints influence on Zr .
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
the length is larger than a fixed number. It should be noticed that constraints’ influence on SNR is quite different under different noise situations, such as in white noise case, the baseline constraint degrades the SNR more badly than in this case. 5. Conclusion The CDMF method is an effective way to calculate the optimum filter with arbitrary noise and constraints. It has been shown that the arbitrary precise flat-top width can be easily get and baseline disturbances with frequency from dc to about several hundreds Hz can be well removed. In company with pulse shape estimator by pole-zero identification method [9], or by least square fitting [12], the CDMF method is very suitable for adaptive, self-calibrating digital spectroscopy which is user friendly and flexible to varying operative conditions.
351
References [1] S. Riboldi, R. Abbiati, A. Geraci, E. Gatti, IEEE Trans. Nucl. Sci.NS- 52 (4) (2005) 954. [2] S. Buzzetti, M. Capou, C. Guazzoni, A. Longoni, R. Mariani, S. Moser, IEEE Trans. Nucl. Sci.NS- 52 (4) (2005) 854. [3] M. Bertolaccini, C. Bussolati, E. Gatti, Nucl. Instr. and Meth. 41 (1966) 173. [4] V. Radeka, N. Karlovac, Nucl. Instr. and Meth. 52 (1967) 86. [5] M.O. Deighton, IEEE Trans. Nucl. Sci. (1969) 68. [6] A. Pullia, G. Ripamonti, Nucl. Instr. and Meth. A 391 (1997) 301. [7] E. Gatti, A. Geraci, G. Ripamonti, Nucl. Instr. and Meth. A 381 (1996) 117. [8] E. Gatti, A. Geraci, S. Riboldi, G. Ripamonti, Nucl. Instr. and Meth. A 523 (2004) 167. [9] A. Pullia, IEEE Trans. Nucl. Sci.NS- 48 (4) (2001) 1234. [10] G. Strang, Linear Algebra and Its Application, second ed., Academic Press, New York, 1984. [11] A. Geraci, M. Zambusi, G. Ripamonti, IEEE Trans. Nucl. Sci.NS- 43 (2) (1996) 731. [12] G. Ripamonti, A. Geraci, Nucl. Instr. and Meth. A 400 (1997) 447.
Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351 www.elsevier.com/locate/nima
Constrained digital matched filter method for optimum filter synthesis Wen Xiangyang, Wei Yixiang Department of Engineering Physics, Tsinghua University, 212, Beijing, China Received 29 November 2005; received in revised form 19 December 2005; accepted 19 December 2005 Available online 17 January 2006
Abstract We present a new method to directly calculate the optimum filter in presence of any additive stationary noise, with arbitrary time and domain constraints (flat-top, zero-area, etc.). A more concise re-deduction of digital penalized LMS method (DPLMS) is given. This method is fully developed, and synthesis results of a typical situation are given and compared with the DPLMS method. Optimum filter can be synthesized without a prior knowledge of the noise power spectral density, which makes it suitable to be used in adaptive, selfcalibrating digital spectroscopy. r 2006 Elsevier B.V. All rights reserved. PACS: 29.30.Kv; 02.10.Sp Keywords: Digital filter; Digital spectroscopy; Optimum filter synthesis
1. Introduction Digital pulse processing systems have been broadly used in high-resolution spectroscopy. One of the major causes of inaccuracy in such systems is the electronic noise introduced by the setup, and the signal-to-noise ratio (SNR) can be greatly improved by appropriate digital processing [1,2]. So digital filter plays a crucial role on the system performances, therefore how to synthesize the optimum filter is a key problem in the design of digital spectroscopy systems. When referring to synthesizing ‘‘optimum filter’’, firstly, we should answer the following question: what type does the optimum filter belong to in the sense of SNR? Gatti and Redeka had proved that time-invariant linear filter can achieve the maximum SNR even without the constraint of linearity or time-invariant, in other words, if we find a filter which achieves the maximum SNR in the scope of timeinvariant and linear filters, it is the optimum filter in general [3,4]. It is impossible, at least in principle, to find a better filter even in a more broad scope based on the SNR criterion. Thus in the following we find the optimum filter Corresponding author. Tel.: +86 10 6279 4799.
E-mail address: [email protected] (W. Xiangyang). 0168-9002/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2005.12.199
from time-invariant linear filters. In many situations arising in practice, the optimum filter not only yields the best SNR, but also should satisfies additional constraints imposed from time domain or frequency domain, such as, flat-top needed for minimizing ballistic deficit effect, overall zero-area of the pulse response for baseline rejection [5,6]. There are two types of linear digital filter namely finite impulse response (FIR) filter and infinite impulse response filter. FIR filter uses only current and past input samples to obtain current output value, and this characteristic is very important, since it means that the input pulses occurring N T times ago will have no contribute to the current output, where N is the length of FIR filter and T is the sample interval. Also FIR filter is not sensitive to filter coefficient quantization and easy to implementation. So FIR filters are very suited for digital spectroscopy applications, although they need more computing power, which is not a problem for current DSP processors and FPGA. Extensive researches have been carried out in the past years, and many methods are proposed to synthesize continuous-time weighting functions which cannot be used directly, and digital filters are derived by sampling its analogue counterpart which are sub-optimum [7]. On the
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
other hand, direct synthesis digital filters is more straightforward and easier to deal with. The DPLMS method introduced by Gatti et al. [8] is a general procedure to directly find the optimum digital filters, which implements a relaxed constraints fulfillment technique which allows achieving constraint accomplishments up to the required level of precision, to avoid further and needless degradation in measurement resolution. Many methods are summarized and fully critically reviewed in Ref. [8]. If several constraints exist at the same time, filter synthesized by DPLMS cannot perfectly satisfy all the constraints. In some practical applications, that every constraint is perfectly fulfilled is more preferable, we propose a new method to meet this need, and how constraints are imposed on the filter is explained in detail. 2. Constrained digital matched filter method
347
point N is, Ps ðNÞ ¼ jhT sj2 .
(2)
Since noise is stationary, output noise power is the same at every point, Pv ðNÞ ¼ EjhT vðN : 1Þj2 ¼ hT EfvðN : 1ÞvðN : 1ÞT gh ¼ hT Rv h
ð3Þ
where Rv is the noise autocorrelation matrix, 0 1 rð0Þ . . . rðN 1Þ B C .. .. .. C Rv ¼ B . . . @ A rðN 1Þ rð0Þ where rðnÞ is autocorrelation function, Rv is positive definite. The SNR at point N is jhT sj2 hT R v h
Constrained digital matched filter (CDMF) method is based on the idea that the optimum FIR filter should maximize the SNR function, and at the same time, satisfy all the constrains imposed on by linear equations. Thus the search of the optimum filter is equivalent to solving an equality-constrained optimization problem.
where SNRðhÞ is a function of h.
2.1. Derivation of SNR
A FIR filter’s output is the convolution of input signal with its impulse response, the output is,
The input sequence obtained by sampling the output of the preamplifier can be written as,
ysðnÞ ¼ sðnÞ hðnÞ ¼
xðnÞ ¼ sðnÞ þ vðnÞ where sðnÞ is the input signal sequence with known shape, and vðnÞ is input zero-means noise sequence which is uncorrelated with sðnÞ. The input signal shape can be get as a prior knowledge, or by sampling and averaging [8] or by pole-zero identification method [9]. The CDMF method uses only noise autocorrelation function, which can be calculated from the power spectrum and system functions before the ADC, or as in many situations arising in practice that we have no exact knowledge about the noise, can be estimated with Eq. (1) by sampling a long noise sequence. rðmÞ ¼
1 N jmj
N1m X
xðnÞxðn þ mÞ.
(1)
n¼0
The output of the N-taps FIR with impulse response h ¼ ðhð0Þ; . . . ; hðN 1ÞÞ is, yðnÞ ¼ hðnÞ xðnÞ ¼ hT sðn : n N þ 1Þ þ hT vðn : n N þ 1Þ. From the input signal sðnÞ we choose a N-point segment which contains the maximum power among all the N-point long segments. For convenience we suppose that this segment starts from point one, and define s as s ¼ ðsðNÞ; . . . ; sð1ÞÞ. Among the output sequence, point N should have the maximum SNR. The signal power at
SNRðhÞ ¼
(4)
2.2. Time domain constrains
N 1 X
hðkÞsðn kÞ.
(5)
k¼0
From Eq. (5) time domain constraints on the output signal shape can be easily applied. Here we give an example on how flat-top constrain is implemented. If a M point flat-top is required, we force the successive (M+1) output to be equal, that is to say, ysðNÞ ¼ ¼ ysðN MÞ. Thus we get M equations, ysðN M þ 1Þ ¼ ysðN MÞ ysðNÞ ¼ ysðN MÞ. Substitute ysðnÞ with Eq. (5), and rewrite it in matrix form, 0 1 xðN m þ 1Þ xðN mÞ . . . xðm þ 2Þ xðm þ 1Þ B C .. .. .. B C B C . . . @ A xðNÞ xðN mÞ xð1Þ xðm þ 1Þ 1 0 hð0Þ 0 1 0 C B C B .. C B C B . C ¼ B ... C B C @ A B B hðN 1Þ C A @ 0 namely, At h ¼ bt .
(6)
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
348
So the flat-top constraint is converted to constraint of linear equations about h. 2.3. Frequency domain constrains
with the null space of AðLT Þ1 , according to the relation between null space and range space, H is the orthogonal complement of range of L1 AT [10]. We performs the QR decomposition on L1 AT , L1 AT ¼ QR ¼ ðq1 ; . . . ; qm ; qmþ1 ; . . . ; qN Þ R. |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl}
The frequency response of a FIR filter is,
T
HðoÞ ¼ DTFTðhðnÞÞ ¼
N1 X
RðA L Þ
hðkÞe
jko
(7)
k¼0
which is fully decided by hðnÞ, thus any frequency domain constraints on the filter are transformed to constraints on the filter coefficients. Here we give an example on how impose constraint on the filter to remove baseline shift. The meaning of removing baseline is that output will be zero if input stay on a constant level, that is to say, the frequency response at 0 Hz is 0, from Eq. (7) we get, N 1 X
hðkÞ ¼ 0
k¼0
rewrite it in matrix form, Af h ¼ bf
(8)
where Af ¼ ð1; . . . ; 1Þ; bf ¼ 0. We combine the foregoing two constraints as follows: Ah ¼ b 0
(9) At
1
B C A ¼ @ A f A;
0
bt
1
The columns of Q can be divided into two groups: the first M columns of Q form an orthonormal basis of RðAT L1 Þ, since RðAT L1 Þ is the orthogonal complement of subspace H, so the last N M columns of Q form an orthonormal basis of H, H ¼ ðqmþ1 ; . . . ; qN Þ. Thus we get the projection matrix onto H, PH ¼ HH T .
(15) L1 v s
onto H, the residual The project matrix PH projects part is orthogonal with any vector in H, thus L1 v s is divided into two parts: one equal to PH L1 s lying in H, the v other is orthogonal with H, therefore is orthogonal with LTv h. We get, T T 1 ðLTv hÞT ðL1 v sÞ ¼ ðLv hÞ ðPH ðLv sÞÞ.
Eq. (11) can be rewritten as, SNRðhÞ ¼
2 jðLTv hÞT PH ðL1 v sÞj ðLTv hÞT ðLTv hÞ
(16)
which, according to the Cauchy–Schwartz inequality, attains its maximum
B C b ¼ @ bf A ¼ 0
2 SNRmax ¼ kPH L1 v sk
when the filter satisfies
where A has full row rank.
(17) LTv h
¼
cPH L1 v s,
h ¼ cðLTv Þ1 PH L1 v s,
2.4. Solve the CDMF problem From Eqs. (4) and (9), we get an equality-constrained maximization problem, namely, hT ssT h max SNRðnÞ ¼ T h2RN h Rv h
(14)
H
1
or equivalently (18)
where c is any constant. The optimum filter can be scaled in any desirable way. 2.5. Remarks on CDMF method
s:t:
Ah ¼ 0.
(10)
which, according to the Cauchy–Schwartz inequality, attains its maximum when the filter satisfies
If the noise is white, Eq. (13) becomes to h ¼ cs, the filter is a scaled time-reverse replica of the input signal shape, this property resulted in the term matched filter. Eq. (13) provides the optimum filter for color additive noise. The optimum filter with constraints is synthesized by three cascaded steps: whitening, projection and matching: whitening step transforms the noise autocorrelation Rv into a unit matrix, and projection step projects the whitened input signal onto the space determined by the constraint equation, and matching is the last step.
LTv h ¼ cL1 v s
3. New deduction of DPLMS method
Perform Cholesky decomposition on Rv , we get, Rv ¼ Lv LTv . Rewrite Eq. (4) as, SNRðhÞ ¼
2 2 jhT Lv L1 jðLTv hÞT ðL1 v sj v sÞj ¼ T T h Lv LTv h ðLTv hÞ ðLTv hÞ
(11)
(12)
where c is any constant. Without the constraints, we get 1 h ¼ cðLTv Þ1 L1 v s ¼ cRv s.
(13)
We define H as the subspace spanned by all possible LTv h, since h is constrained by Eq. (9), L1 v s may not lie in H, in other words, Eq. (12) may have no solution. H is identical
To begin with, the main idea of this method is recalled: the optimum filter is found by minimizing the weighted sum of the output noise variance and other terms that penalize constraint violations, all penalty terms are in quadratic form. In case of only flat-top and baseline
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
rejection constraints are concerned, the problem can be expressed as, 8 T 8 PV ðNÞ ) 0 > h Rv h ) 0 > > > > > T > < < Y ðNÞ ) 1 kh s 1k2 ) 1 ¼) > > A t h ) bt kAt h bt k2 ) 0 > > > > > : : kA h b k2 ) 0: A f h ) bf f f Assume the weights are ½a1 ; a2 ; a3 ; a4 , the objective function can be written as, FðhÞmin ¼ a1 hT Rv h þ a2 khT s 1k2 þ a3 kAt h bt k2 þ a4 kAf h bf k2 ¼ hT Ah 2bT h þ a2 þ a3 bTt bt þ a4 bTf bf T
a3 ATt At
a4 ATf Af ,
ð19Þ T
where A ¼ a1 Rv þ a2 ss þ þ b ¼ a2 s T þ T T a3 bt At þ a4 bf Af . Since Rv is positive definite and symmetric and other terms are real and symmetric, so A is nonsingular, the solution of Eq. (19) is, DFðhÞ ¼ Ah b ¼ 0 ) h ¼ A1 b.
(20)
349
eliminates the need for any flat-top, as is the case of the present example. A special point should be noticed is that, in digital systems, a flat-top further by one sampling time is needed because of the random jitter between sampling comb and pulse arrival time [11]. In the following, a 15 points flat-top constraint is assumed which is not really needed for the present detector, but just to exhibit the capability of flat-top filters synthesizing. The optimum filters for four cases are synthesized respectively, A: no constraint, B: flat-top constraint, C: baseline rejection constraint, D: flat-top and baseline rejection constraint. The synthesis results are plotted in Fig. 1, it can be seen that output A in Fig. 1(a) is more cusp-like than its counterpart in Fig. 1(b), which is in accord with the matched filter theory. 4.3. Comparisons between DPLMS and CDMF It is easy to verify from Eq. (18) that the synthesis result strictly satisfies all the constraints, so the CDMF method
4
x 104 A B C D
4. Synthesis results and discussions 3 2 Output
In the following, the CDMF method described in Section 2 was applied to a typical situation, the DPLMS method re-deduced in Section 3 was also applied for comparison.
1
4.1. Input conditions
0
We take the small silicon detector as an example, this simulated system has the following conditions: d input current pulse, only series noise and parallel noise with noise corner time tc . The sample frequency is 10 MHz, and the analog system before the ADC is the cascade of three parts: charge-sensitive amplifier, whitening differentiator with time constant followed by a simple R-C anti-aliasing filter, the transfer function is, K ðs þ 1=tc Þðs þ 1=ta Þ
pa2 jtj=ta e 2ta
0
10
20
30
(a)
40 50 time(n)
60
70
80
90
8 A B C D
6
(21) 4
where tc ¼ 0:5 ms; ta ¼ 1:0 ms, tc and ta is the time constant of the differentiator and the anti-aliasing respectively. The overall autocorrelation function is, RðtÞ ¼
−2
(22)
where a2 ¼ 1:52e015 v2 s.
2 Output
F ðsÞ ¼
−1
0 −2 −4
4.2. Synthesis results −6
When designing optimum filter for a given system, constraints should be imposed according to the specific conditions, in cases where the charge collection time of the detector can be regarded to as negligible, which in practice
(b)
0
10
20
30
40
50
60
70
80
90
time(n)
Fig. 1. Comparison of the weight functions with different noise situations: (a) color noise described by Eq. (22), (b) white noise.
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
350
can get perfect flat-top and absolutely baseline rejection at the same time. While from Eq. (20), if no constraints exist, the DPLMS method leads to matched filter and gets the same result as CDMF. If one or more constraints are imposed on, the result is a compromise between the SNR and fulfillment of the constraints. If there is only one constraint, by emphasizing the associated weight coefficient the constraint can be almost perfectly satisfied, at the expense of a degraded SNR, which can be seen from Fig. 2(a) while the flat-top weight value changes from less effective ða3 ¼ 30Þ to highly effective ða3 ¼ 1000Þ, the top of the output is more and more flat. However, if there are more than one constraints, the constraints cannot be perfectly satisfied at the same time, which can be seen from Fig. 2(b), the input sequence staying on a dc voltage of 2, is filtered by filters synthesized with different a4, the output shift is 0.2 when a4 ¼ 1, and is almost completely removed when a4 ¼ 10, but while the ability of removing baseline enhanced, the flatness of the output top is worsened.
1.0015
α3 = 30 α3 = 100 α3 = 1000 CDMF
1.001
1.0005
1
Theoretic SNR for several situations is computed for the filters synthesized by the two methods, with the foregoing input conditions, and is normalized to the SNR given by CDMF with no constraint. Some result is given in Table 1. When both the constraints exist, the two methods give almost the same SNR, but, in practice applications, the imperfect fulfillment of the constraints by DPLMS would worsen the SNR, which is also proved by simulation. 4.4. Constraints’ influence on SNR We define the absolute amplitude-to-rms noise ratio as Z ¼ SNR1=2 , from Eq. (17), the possible maximum Z with given input s is, Z1 ¼ limN!1 kL1 v sk, for convenience, we define the relative amplitude-to-rms noise ratio as Zr ¼ Z=Z1 . When synthesizing the optimum filter, the first decision we have to make is how long the FIR filter is, for this purpose we calculate Zr as a function of FIR length N, the result shows that relatively short FIR length is sufficient to achieve SNR fairly close to the ideal SNR Z1 (Zr 40:99 with FIR length 30). Constraints imposed from the time domain or frequency domain reduce the filter freedom degrees, therefore would results in a degraded SNR. The influence on Zr of various filter length and constraints is plotted in Fig. 3. There are six curves in this figure. In fact, because the three curves of filter length 30; 50; 80, respectively, without baseline rejection are almost identical, which are labeled as ‘Non-BLR’ as a whole, we can conclude that without baseline rejection constraint, filter length has little impact on the SNR when Table 1 The SNR of DPLMS with both constraints is 0.626, the SNR of DPLMS is computed for four cases, ½a1 ; a2 is fixed to ½1; 1
0.9995
2
3
4
5
6
7
8
9
½a3 ; a4
½0; 1
½1; 0
½1; 1
½10; 10
SNR
0.876
0.818
0.626
0.626
n
(a) 5
α4 = 1 α4 = 10 α4 = 1000 CDMF
0.4
4
0.2
3
0
1 Non−BLR 30+BLR 50+BLR 80+BLR
0.9
2
0.8
1 0.7
−1
4.6
−2
0.5
4.2 4
0.4
3.8
−4 −5
0.6
4.4
−3
(b)
ηr
0
0.3 0
20
40
60
80
100
120
Fig. 2. Filters synthesized for the following conditions: flat-top width 6, FIR length 60, weight value is 1 if not explicitly asserted: (a) flat-top synthesized with different value of a3, (b) baseline remove effects with different value of a4.
0.2
0
1
2
3
4
5
6 7 8 9 10 11 12 13 14 15 Flat top Length
Fig. 3. Filter length and constraints influence on Zr .
ARTICLE IN PRESS W. Xiangyang, W. Yixiang / Nuclear Instruments and Methods in Physics Research A 560 (2006) 346–351
the length is larger than a fixed number. It should be noticed that constraints’ influence on SNR is quite different under different noise situations, such as in white noise case, the baseline constraint degrades the SNR more badly than in this case. 5. Conclusion The CDMF method is an effective way to calculate the optimum filter with arbitrary noise and constraints. It has been shown that the arbitrary precise flat-top width can be easily get and baseline disturbances with frequency from dc to about several hundreds Hz can be well removed. In company with pulse shape estimator by pole-zero identification method [9], or by least square fitting [12], the CDMF method is very suitable for adaptive, self-calibrating digital spectroscopy which is user friendly and flexible to varying operative conditions.
351
References [1] S. Riboldi, R. Abbiati, A. Geraci, E. Gatti, IEEE Trans. Nucl. Sci.NS- 52 (4) (2005) 954. [2] S. Buzzetti, M. Capou, C. Guazzoni, A. Longoni, R. Mariani, S. Moser, IEEE Trans. Nucl. Sci.NS- 52 (4) (2005) 854. [3] M. Bertolaccini, C. Bussolati, E. Gatti, Nucl. Instr. and Meth. 41 (1966) 173. [4] V. Radeka, N. Karlovac, Nucl. Instr. and Meth. 52 (1967) 86. [5] M.O. Deighton, IEEE Trans. Nucl. Sci. (1969) 68. [6] A. Pullia, G. Ripamonti, Nucl. Instr. and Meth. A 391 (1997) 301. [7] E. Gatti, A. Geraci, G. Ripamonti, Nucl. Instr. and Meth. A 381 (1996) 117. [8] E. Gatti, A. Geraci, S. Riboldi, G. Ripamonti, Nucl. Instr. and Meth. A 523 (2004) 167. [9] A. Pullia, IEEE Trans. Nucl. Sci.NS- 48 (4) (2001) 1234. [10] G. Strang, Linear Algebra and Its Application, second ed., Academic Press, New York, 1984. [11] A. Geraci, M. Zambusi, G. Ripamonti, IEEE Trans. Nucl. Sci.NS- 43 (2) (1996) 731. [12] G. Ripamonti, A. Geraci, Nucl. Instr. and Meth. A 400 (1997) 447.