1760
ANALYTICAL CHEMISTRY, VOL. 51, NO. 11, SEPTEMBER 1979
Efficient Computation of Polynomial Smoothing Digital Filters M. U. A. Bromba and Horst Ziegler” Fachbereich Naturwissenschaften I, Angewandfe Physik, Gesamthochschule Paderborn, Warburger Sfr. 100, D-4790 Paderborn, West Germany
Smoothing spectra with digital filters has many advantages over conventional RC filters. The polynomial smoothing digital filters are especially suitable for spectrometry. This paper presents a method for reducing the computational requirements by several orders of magnitude by using a recursive algorithm. The method is stable and has zero numerical error while preserving all other properties of Its well known nonrecursive ancestors. Formulas are given for the zero, second, and fourth degree polynomial filters and computational and storage requirements are discussed.
Smoothing digital filters are a valuable tool for improving the signal-to-noise ratio of spectrometric measurements. Especially suitable to all types of spectrometry are the finite polynomial smoothing filters, described by Savitzky and Golay ( I ) and others (2-7). They conserve area, symmetry, line position (for symmetric lines), and higher moments while yielding a signal] noise improvement which can be near the theoretical limit of the matched filter. Realtime filtering a t a rate suitable for spectrometry, however, was previously not feasible because of the excessive number of multiplications required and the short time available between samples. This paper presents a new algorithm for implementing these filters. It dramatically reduces the computational requirements to the point that a standard microprocessor can easily implement such a filter in realtime. Theory-Nonrecursive. Starting with the general equation for a digital filter, we have m
D f [ k ]=
C a [ n ] . f [ k- n ] = ( a * f ) [ k ] a=
(1)
-m
where D denotes the (linear) filter operator, f [ h ]the hth signal point (-m < k < m ) , a[n]the nth point of the filter function and “*” the convolution operation, which is both associative and commutative. The polynomial smoothing filter operations of degree 2M conserve all moments of the signal up to the (2M 1)th moment. Examples are: (1) the running average ( M = 0) which conserves area (zero moment) and center of gravity (first moment) and (2) the well known second degree filter which conserves in addition the second moment of a line, which is the true measure of the spectrometric line width. The coefficients a z M [ nfor ] these finite (azhf[n]=0 for In1 > N) filter functions are:
+
dz = (2N - 1)(2N + 1)(2N + 3)/3 d4 = 4(2N - 3)(2N - 1)(2N
(6)
+ 1)(2N+ 3)(2N+ 5)/15 (7)
For precomputed coefficients the number of multiplications required is (N I),N being the half width of these finite filter functions. For achieving high accuracy in spectrometric applications any distorting RC-prefiltering must be reduced to an absolute minimum. This in turn requires a fairly high time resolution in the digital filter (i.e., dense samples) and hence fairly high values for the filter width N . Suggested values for spectrometric applications are a time resolution of 20 ms and a filter width of N = 100 to 500 for a scan speed of 10 s/line. During each 20-ms interval we would have to perform up to 500 multiplications, a task clearly beyond the scope of current microprocessors. Theory-Recursive. Let us introduce the linear difference operator A with the property
+
WI = f [ k + 11 - f [ k l
(8)
If we apply this operator (2M + 1)times on both sides of the filter Equation 1, we can write
a2M+1 (DSMf) =
(9)
(asM*f)
Since the difference operator A is also a convolution operator, thus having the properties of associativity and commutativity, we can rewrite the right hand side of Equation 9 We thus have to compute the (2M + 1)th difference of the original coefficients uzM and convolute this result with the input signal f. Since the filter coefficients uzM form a polynomial of degree 2M, the (2M + 1)th difference will be zero for almost all points. Only near the end points ( + N and -N) there is a contribution because of the “discontinuous” jump to zero of the finite filter function am The number of these nonzero terms which determine the number of multiplications required does not depend on the width N of the filter. In Table I an example of this procedure is given for 2M = 2 and N = 7 (the well known 15-point quadratic smoothing). All values given are without the common normalizing integer denominator d z = 1105. Obviously the convolution of the input signal f with the new coefficients A3az considered in Equation 10 needs less calculation than the original convolution uzM*f in Equation l , since A%, has only 6 nonzero coefficients in contrast to 15 nonzero coefficients of u2. With the coefficients A3uz[n] of Table I, Equation 10 becomes m
A3(D2f,[k]= 1 -(-78f[k 1105
The normalizing integer denominators are:
d, = 2N-t 1
(5)
0003-2700/79/0351-1760$01 .OO/O
n=-m
(A3az)[n]f[k - n] =
+ 101 + 221f[k + 91 - 153f[k + 81 + 153f[k - 51 - 221f [ k - 61
1979 American Chemical Society
+ 78f[h
-
71) (11)
ANALYTICAL CHEMISTRY, VOL. 51, NO. 11, SEPTEMBER 1979
Since we are not interested in the values A3(D2f)[k]but in D . j [ k ] ,we take the third difference of D2f and get
A3(Dzf,[h]= D j [ k + 31 - 3Dd[h + 21
+ 3D$[h + 11
Inserting this into Equation 11, replacing k by k solving for D d [ k ] ,we obtain
-
3, and
D j [ k ] = D.j[h - 31 - 3 ( D j [ h- 21 - D j [ k - 11) 78 -(f[h + 71 - f [ h- 101) + -221 (f[k + 61 - f [ k 1105
1105
91)
153 -(f[h+ 1105
-
51 - f [ h - 81) (13)
Equation 13 may be viewed as a recursive digital filtering process (see (6)) with finite memory and having exactly the same effect on a spectrometric signal as the nonrecursive 15-point quadratic Savitzky-Golay smoothing. Now let us generalize the example for an arbitrary width N of the filter and an arbitrary degree 2M. We apply the binomial expansion to the left hand side of Equation 10 and rewrite the right hand side noting which terms of the convolution sum will be nonzero. We then replace h by k - (2M + 1) as we have done in the example and solve for D2.vf[k]. 2M
D&[h] =
(2Mm+ m=O
I)
(-1)2M-"D&[h - 2M - 1 + m] +
2M
2 AZM+la2M[N m]Cf[h -
-
N
-
2M
-
1
m=O
-
m ] ) (14)
A smoothed point a t position k ( D P M f [ kthus ] ) depends both on the previous filtered values (first sum) and a few unfiltered values (second sum). The resulting formulas for the most common filters are: -
fth
-
N
-
11)/dO (15)
D d [ h ]= D j [ h - 31 - 3 ( D j [ k - 21 - D j [ h - 11) (2N - 1 ) ( N - l ) ( f [ k+ N ] - f [ k - N - 3 ] ) / d , + ( 2 N + 3 ) ( 2 N - 1) ( f [ h + N - 11 - f [ h - N - 2 ] ) / d , ( 2 N + 3 ) ( N + 2 ) ( f [ h+ A - 21 - f [ k - N - l ] ) / d , (16) D j [ k ] = D j [ h - 51 - 5 ( D j [ h- 41 - D j [ k - 11) + 1 0 ( D J [ h - 31 - D$[h - 21) (2N - l)(2N - 3 ) ( N 2)(N - l ) ( f [ h N] - f [ h - N - 5 ] ) . 2 / d , - (2iV + 5 ) ( 2 N - 1)(2N- 3 ) ( 2 N -4) X ( f [ h + N - I ] - f [ h - I L T ~4 ] ) * 2 / d ,t ( 2 N + 5 ) ( 2 N 3)(6W+ 6N - 2 2 ) ( f [ h+ N - 21 - f [ h - N - 3]).2/d* - (2N + 6 ) ( 2 N + 5 ) ( 2 N + 3)(2N - 3 ) X ( f [ h + N - 31 - f [ k - N - 2]).2/d4 (2N + 5 ) ( 2 N 3 ) ( N + 3 ) ( N + 2 ) ( f [ h+ N - 41 - f [ h - N - l ] ) * 2 / d 4 (17)
+ +
+
1105. n
+
Obviously the number of multiplications now only depends on the degree 2M and not on the width N. Numerical Stability. Using recursive techniques for finite filters usually leads to a numerical instability of the algorithm. Since either the coefficients themselves or the products have t o be rounded to a finite word length, a numerical error is introduced. Because of the recursive nature of the formula, this error propagates through all following computations, possibly leading to a constant offset or even to a divergence of the algorithm. If we analyze the recursive filter equations, 15 to 17, we note that the coefficients for the unfiltered input values are rational numbers which cannot be represented
[ n1
QZ
- 11
1105.
1105,
- 78
143
-4
-78 -13 42 87
65 55 45 35
-3
122
25
-2
147
-1 0 1
162
13 - 78
15 5 -5 - 15 - 25 - 35 -45 - 55 - 65 78
- 10 - 10 - 10 - 10 - 10 - 10 - 10 - 10 - 10 - 10 - 10 - 10 - 10
0 0
0 0
-7 -6 -5
2
167 162 147
3
122
4 5 6
42
7
0 0 0
87 -
8
9
1105,
( A a , ) [ n1 (A'Qz)[nI l:A3Q,) [ n1
0 0 0 0
- 10 -9 -8
0 0 - 78
143 - 78
0 0
0
78 221 - 153 -
0 0 0 0 0 0 0 0 0 0 0 0
153 - 221 78 0 0
Table 11. Comparison of Computational and Memory Requirements for Recursive/Nonrecursive Polynomial Filters
+ m] - f [ h + N
DOf[hl = Doftk - 11 + ( f [ h + NI
Table I. Coefficients and Differences between Coefficients for a 15-Point Quadratic Smooth
-
D d k I (12)
1761
nonrecursive memory data points coefficients computation mu1tiplications addlsu b trac t
recursive
2N+ 1 N L 1
:!N :IM
N+ 1 21v 4- 1
:IM tIM
+ 4M + 3 +
1
+
1
f
2
without error in fixed word length. We note, however, that they all share a common integer denominator. We now multiply both sides of the equation with this denominator dZM and define E2Mf[kl
= d2MD2Mf[kl
(18)
We then store these multiplied partial sums E 2 v f [ k ]instead of the original D2,f[k]. Thus all coefficients, all data, and all intermediate results are pure integers, which can be represented error-free in a finite word length. To output the final filtered value, we divide the last recursive sum E z ~ f [ k ] by the denominator dZM. The result of this normalizing division, however, is only used to feed the D/A converter and is not used in any recursive calculation. To be sure that the start of the recursion process does not introduce any instabilities, one must carefully choose a set of initial conditions. This may be done by setting all values f [ k ] ,h < ko, preceding the first spectrometric signal value f [ k o ] as well as DzMf[kO - N - (2M + 111. . . DLvf[kO- N - 11 equal to a constant (e.g., zero or f [ k o ] before ) starting the calculation of DZMf[kO - Nl. This leads to a stable recursive finite filter with zero computational error. Comparison. Table I1 compares the number of multiplications and additions/subtractions for the two methods. Also the memory requirements for both data points and coefficients are given. I t is assumed that pairs of data with the same multiplier are first added,/subtractchd and then multiplied. When comparing computational requirements, we must also take into account the required word length of the multiplications involved. Comparing the biggest coefficients in both the nonrecursive and in the recursive formulation (Equations
1762
ANALYTICAL CHEMISTRY, VOL. 51, NO. 11, SEPTEMBER 1979
2 to 4, respectively, 15 to 17) we note that they are both of order p. This corresponds to the figures of Table I giving a biggest coefficient of 167 for the nonrecursive vs. 221 for the recursive case. For a typical situation with N = 500 and 2M = 2 (quadratic filter), the critical number of multiplications drops from 501 to just 4, a dramatic saving which even exceeds that of the famous Fast Fourier Transform algorithm. Thus, except for extremely narrow filter functions (N < 3M), it is always better to use the recursive formulas. Implementation. While the normal nonrecursive smoothing algorithm can be easily implemented using real numbers in a FORTRAN program, the stability requirements of the recursive filter demand integer arithmetic for all calculations except the final normalizing division. This cannot be done effectively with a FORTRAN program where integer word length is usually very limited. Any standard microprocessor has the required capability
for multibyte integer arithmetic. We have successfully implemented such a realtime digital filter based on this new recursive technique using a standard 8085A-microprocessor. It can handle all these filters for up to 2M = 4 (quartic filter) and any filter width up to N = 1800 with a realtime sampling interval of only 20 ms and an accuracy which is limited only by the AJD and D/A conversion accuracy.
LITERATURE CITED (1) Savitzky, A.; Golay, M. J. E. Anal. Chem. 1964, 36, 1627-1638. (2) Steinier, J.; Termonia, Y.; DeRour, J. Anal. Chern. 1972, 4 4 , 1906-1909. (3) Madden, H. H. Anal. Chem. 1978, 50, 1383-1386. (4) Ernst, R . R . Adv. Magn. Reson. 1966, 2, 1-135. (5) Wiilson, P. D.; Edwards, T. H. Appl. Speclrosc. Rev. 1976, 72,1-81. (6) Hamming, R. W. "Digital Filters"; Prentice-Hall: Englewood Cliffs, N.J., 1977. (7) Bromba, M. U. A., Diplomarbeit; Paderborn, 1978. (8) Blum, M . IRE Trans. Inf. Theory 1957, IT-3, 178-182.
RECEIVED for review November 20, 1978. Accepted May 7 , 1979.
Trace Analysis of Lanthanides by Laser Excitation of Precipitates Frederick J. Gustafson and John C. Wright" Department of Chemisfry, University of Wisconsin, Madison, Wisconsin 53706
Selective excitation of probe ion luminescence (SEPIL) has been applied to the trace determination of nine lanthanides without prior chemical separation. The lanthanides are quantltatively coprecipitated in calcium fluoride and their 4f electron solid-state luminescence is excited by a tunable dye laser. The selectivity, detection limlts, concentration dependence, extent of coprecipitation, sample preparation parameters, and the temperature dependence of the emission intensity of the method have been investigated. The interference effects of 36 anions and cations are presented.
In a previous publication, a new analytical technique was introduced for trace analysis by selective excitation of probe ion luminescence (SEPIL) ( I ) , in which intense fluorescence was excited from lanthanide ions coprecipitated in calcium fluoride by a pulsed, narrow bandwidth, tunable dye laser. Selective excitation of different crystallographic sites of the same ion is possible because the lanthanides have narrow energy levels when doped into ordered crystalline media (see Reference 1 for a brief description of the theory of lanthanide solid state luminescence). The lanthanide 4fn electron energy level structure is different for each particular lanthanide ion making it possible to excite each lanthanide selectively in the presence of other lanthanides. This work will report the extension of SEPIL to the determination of nine lanthanides, praseodymium, neodymium, samarium, europium, terbium, dysprosium, holmium, erbium, and thulium in the same sample. All of these ions fluoresce and quantitatively coprecipitate in calcium fluoride. A single crystallographic site per ion (labeled the G1 site) is found to dominate the spectra of all nine lanthanides under proper 0003-2700/79/035 1- 1762$01 0010
ignition conditions of the calcium fluoride. The peak height intensities of the G1 site for each lanthanide is proportional to the original concentration of the lanthanide in solution before precipitation. The detection limits for the lanthanides were experimentally determined and range from 0.4 pg/mL of Eu3+to 450 pg/mL for Nd3+. The detection limits for most of the lanthanides are lower than published detection limits for any other analytical technique (2-12). This work will also report the results of precipitating the lanthanides in the presence of other ions such as might be encountered in a real sample.
EXPERIMENTAL Apparatus. A complete description of the nitrogen laser pumped dye laser system used to study precipitates was presented by Miller e t al. (13) and will not be repeated here. Reagents. All reagents used in this study were analytical reagent grade and used as supplied. Solutions of 1.2 M Ca(NO& and 0.3 M NH4F were prepared by simply weighing and dissolving the undried salts. The molarities of the calcium and fluoride solutions were shown to be accurate within 0.3% of each other by standardization ( I ) . Approximately 0.03 M solutions of praseodymium, neodymium, samarium, europium, terbium, dysprosium, holmium, erbium, thulium, and ytterbium were prepared by weighing either the oxides or the hydrated chlorides (99.970, Research Organic/Inorganic Chemicals, Belleville, N.J.). The chlorides were dissolved in M "OB to prevent the formation of rare earth hydroxides. The oxides were dissolved in 30 mL of 3 M "OB, evaporated to dryness on a hot plate, and M "OB. All of the lanthanide ion solutions redissolved in were standardized by the EDTA titration method described by Haskin et al. (14). Dilutions of the lanthanide ion solutions were M "OB. Fresh dimade from these stock solutions into lutions were made before all concentration dependence experiments. To minimize the use of unstable low molarity lanthanide ion solutions, small volumes (I mL or less) of concentrated soC 1979 American Chemical Society
