Anal. Chem. 1996, 68, 1054-1057
Anomalies of Convolutional Smoothing and Differentiation K. Y. Hui and M. Gratzl*
Department of Biomedical Engineering, Case Western Reserve University, Cleveland, Ohio 44106
The original Savitzky-Golay convolutional smoothing method for simultaneous data smoothing and differentiation has been a useful tool for processing experimental data for nearly 30 years. However, unlike splining techniques, it leads to truncation of the end sections due to the center-point moving window approach. For a long time this was considered to be a necessary evil that could not be overcome. In 1990, Gorry addressed the truncation problem for uniformly spaced data in a mathematically elegant way. We implemented and tested this algorithm in our laboratory and observed some anomalies in the results produced. We offer explanations for these anomalies and suggestions of modifications so that the observed undesirable effects can be avoided. In 1990 Gorry gave the experimental world a modified Savitzky-Golay1 convolutional smoothing algorithm that would not lead to end region truncation of the experimental data.2 A later augmented version was published to cover nonuniformly spaced data.3 The implementation is simple and straightforward. The authors tested the algorithm for uniformly spaced data on some experimental results and were pleased by its efficiency and elegance. However, there were also some unexpected (at the time) peculiarities produced by the smoothing operation. Upon further investigation, we realized that there are some inherent properties of the least-squares polynomial convolutional smoothing approach that even the new general method did not take into account. We share our findings in the following. The obvious reason that one does not try to fit just one polynomial of a very high order to a whole set of data is that this approach will be computationally inefficient and the results will not even be satisfactory for most cases. For example, polynomials are not well suited to fit curves containing peaks or with asymptotic behavior. At the same time, the piecewise moving window approach allows for a good fit with a relatively low order polynomial. The reason is the distinct adaptive characteristic of this technique which, with proper tuning, can be a very powerful and practical tool for data processing. The intrinsic requirement for the convolutional piecewise polynomial smoothing approach is therefore the dependence on both past and “future” information to obtain a good enough evaluation of the new information. This is why the original Savitzky-Golay1 method truncates at the end sections of the data. What the new method2,3 does is to fit the first 2m + 1 points of a data set on which a 2m + 1 wide window is operating with a (1) Savitzky, A.; Golay, M. J. E. Anal. Chem. 1964, 36, 1627-1639. (2) Gorry, P. A. Anal. Chem. 1990, 62, 570-573. (3) Gorry, P. A. Anal. Chem. 1991, 63, 534-536.
1054 Analytical Chemistry, Vol. 68, No. 6, March 15, 1996
polynomial of the order appropriate for the entire data set. The resulting first m points are then used for the portion that would have been truncated by the Savitzky-Golay method. The last m points are obtained similarly. The algorithm produces an n × n (where n ) 2m + 1) square matrix of convolutional weights. Each of the left-side m columns of weights is used for preserving each of the corresponding leftmost m points of the signal. Likewise, each of the right-side m columns of weights is used for preserving each of the corresponding rightmost m points of the signal. The center column of weights is then used for all the data points in the center. Although this method was arrived at in an elegant way and provided a simple implementation that is computationally very efficient, we have realized that both end section fittings thus resulted are oversimplistic. From a digital signal processing point of view, each column of weight coefficients represents a unique finite impulse response filter which have distinct frequency response characteristics. Therefore, each of the end region points is processed by a unique filter especially suited for that particular point at that particular position in relation to its neighboring points. Of course, if one analyzes each column of weights separately, one would find that they are all low-pass filters with different frequency responses. The center column weights are symmetrical (see tables in ref 2) because the fitting done in the center section is a symmetrical operation: the processed point is always in the center of the digital filter window. Hence, one would expect a linear phase response for the center column of weights. The side columns of weights corresponding to each end section point (again, see tables in ref 2) are all asymmetrical because of the asymmetries of each filter operation: the processed point is never in the center and is never in the same position twice. Hence, one would expect nonlinear phase responses for all these columns of weights. The net consequence of this nontruncating algorithm pertaining to the end sections can be readily identified from empirical results as that each end section is fitted by one single polynomial of the same order chosen for the center section, without the benefit of a sliding window as does the center section. This can be most easily understood by processing an n-point data set using an n-point filter (of size n × n, of course). The end sections thus obtained are more flexible than outright truncation, but they still do not have the flexibility of the center section provided by the adaptivity due to the sliding window. Thus, when the overall filter is tuned for the center section of a raw data set, the polynomial order chosen for this major part of the signal will inevitably be too low for the end sections. In this work we first demonstrated the end section artifacts by using simple simulations. A possible remedy by using higherorder polynomial fitting for the end sections than the center 0003-2700/96/0368-1054$12.00/0
© 1996 American Chemical Society
Figure 1. A simple demonstration of end section truncation. (A) An ideal quadratic curve of y ) x2 (x from -1 to 1 with 0.005 step). (B) Smoothed first derivative of the fundamental curve A using a 75point filter generated by the original nontruncating algorithm (of course, if one has chosen a second-order polynomial, there will not be any truncation). However, the fact that the center portion is correctly processed shows that, for real experimental data (which are rarely ideal polynomials), what is right for the center may be too rigid for the end sections in the original algorithm.
section was tested. Another approach for real-time processing by using only the right-side m columns for fitting the most current m points successively was also tested. EXPERIMENTAL SECTION Three separate experiments were performed. First, a simulated quadratic curve without any noise added was “processed” using a 75-point filter (large filter width chosen to make the end region artifacts clearer) with linear fit and firstorder differentiation (Figure 1). We then determined the end point of a titration curve (Figure 2A) by detecting the zero-crossing of the second derivative of the fundamental data array. For detailed experimental conditions, see refs 4-6. A 25-point filter according to the nontruncating algorithm in ref 2 was first used to compute the smoothed second derivative of an entire titration curve that includes data points well past the end point (curve a in Figure 2B). The same filter was then used to filter the same titration curve up to the end point (curve b in Figure 2B). Then, different orders of polynomial fitting for the center and end sections, third and fifth, respectively, were used for the same titration curve (Figure 2C). In the third experiment, after the initial m + 1 (here, m ) 12) points, each most recent m point after every slide of the window thereafter was computed using only the right-side m columns of the precalculated weight coefficients for a higher-order (fifth in this case) polynomial fitting (Figure 2D). This simulated the realtime processing condition. RESULTS AND DISCUSSION Simulated Data. Figure 1 shows a simulated quadratic curve and its “smoothed” first derivative using linear fitting. Although the first derivative of a linear curve should produce a constant, (4) Gratzl, M. Anal. Chem. 1988, 60, 484-488. (5) Gratzl, M. Anal. Chem. 1988, 60, 2147-2152. (6) Diefes, R. S.; Hui, K. Y.; Dudik, L.; Liu, C. C.; Gratzl, M. Sens. Actuators B, in press.
the center section is correctly computed due to the additional flexibility of the moving window approach, while the end sections are flat. This will always happen when the order of derivative needed is equal to the order of the chosen polynomial. Of course, with the a priori knowledge that we have about the input signal, we should have chosen a second-order polynomial fitting, but real signals are rarely ideal polynomials. The whole objective of filtering is to exclude as much noise as possible without losing useful information. The cutoff frequency of a piecewise polynomial filter correlates with the order of polynomial. Inasmuch as matching the highest signal frequency, one should avoid using higher-order polynomials which will lead to less effective noise filtering. Thus, assuming no a priori information as in realistic cases, the filter in this example is properly tuned for the center section. Unfortunately, this choice is too simplistic for the end sections because of the relative lack of filter adaptivity in the end sections compared to the center portion. Postprocessing Experimental Data. In a conventional leastsquares piecewise polynomial filter (Savitzky-Golay convolutional approach), every point except the leftmost and the rightmost m points of a data array is smoothed with the information of a unique new point in addition to the “old” information, while deleting the oldest point. This results in an adaptive behavior which updates quickly and accurately. The major drawback of this approach is the truncation of the end sections leading to loss of information. To remedy this condition, the nontruncating algorithm2 allows for the computation of additional convolutional weights for each of the end section points. The end section fitting is done with the same set of raw data for every point albeit a different set of weights is used to arrive at each processed point. Therefore, when one properly tunes the filter to smooth the center portion of the data set, one necessarily chooses a polynomial too simple to fit the end sections. For the center portion, one needs to choose a simple enough polynomial that may seem insufficient for the frequency content of the signal because of the adaptive property of a moving (convolutional or spline) filter. This polynomial then turns out to be too rigid for the end sections where only one polynomial curve is fitted to all the raw data points. One of the many important potential applications of a signal processing algorithm without truncation would be in voltammetry. The fact that the modified algorithm overcomes the problem of end section truncations would be significant in processing data without truncating the entire useful voltage window beyond which other interfering processes may start taking place. However, by choosing a polynomial order that enables one to properly preserve the information in the end regions, one might retain excessive noise in the center region, where important information may also reside. The titration curve up to the end point in Figure 2A was used as the input signal to a slightly modified filter. Two different orders of polynomial fitting were used to compute the smoothed second derivative of this curve. The center portion was processed using the center column of weights from an n-point (here, n ) 25) third-order filter and the right-end section was processed using the right-side m (here, m ) 12) columns of weights from another n-point fifth-order filter (Figure 2C), thus resolving the dilemma discussed above (as illustrated in Figure 2B and Figure 2D). Real-Time Signal Processing. Another practical application will be in real-time signal processing such as in the determination of titration end point in a stationary diffusional microtitration Analytical Chemistry, Vol. 68, No. 6, March 15, 1996
1055
Figure 2. An example of applying the original nontruncating algorithm to experimental data. (A) A diffusional microtitration curve of the titration of sodium thiosulfate by sodium triiodide. The end point is obtained by locating the zero-crossing of the second derivative in (B). For experimental details, see refs 4-6. (B) The smoothed second derivative of the fundamental curve A (1-Hz sampling frequency, 25-point filter length, cubic fit): (a) Using data collected well past the end point. (b) Using data up to the end point. The failure of the algorithm to detect the end point (zero-crossing) demonstrates the insufficient flexibility of the polynomial in the end section. (C) Results of the modified algorithm proposed for postprocessing: (a) Smoothed second derivative of the fundamental curve A by using 25-point filters of third-order fit for the center and fifth-order fit for the right-end section. The noise in the center region is reduced, and the end point is correctly detected. (b) The last 25-point fit is shown. (D) The results of a modified algorithm proposed for real-time feature detection. Smoothed second derivative obtained using only the right-side columns of weights. While the end point is now correctly detected, the noise level prior to the end point has increased. This further proves that the same filter parameters could not be used for postprocessing the entire signal without either oversmoothing the end sections or undersmoothing the center section.
scheme. Here, real-time detection of the inflection of the monitored titration curve is crucial. In such a scheme, diffusion is employed as the sole driving force for reagent delivery where each titration has to begin with the same initial conditions regarding the concentration profile of reagent in the diffusion membrane. Therefore, it is important to detect the end point in real time so that the current titration can be terminated immediately and the next determination can be started. Thus, a steady-state reagent concentration profile in the diffusion membrane can be ensured at all times.4,5 The experimental procedures for performing the microtitrations can be found in refs 4-6. Curve b in Figure 2B demonstrates the limitation of the original nontruncating algorithm in this case, which may lead to erroneous interpretations of analytical results. Potential errors demonstrated here call for educated cautiousness. With improved understanding, modifications are possible. Instead of using the whole algorithm as it was originally intended, one can use the different precalculated convolutional weights of 1056
Analytical Chemistry, Vol. 68, No. 6, March 15, 1996
the right-side end section for efficient real-time polynomial fitting. This is like splining except that it is more computationally efficient (hence more suited for real-time processing) because of the convolutional approach. Thus, one should choose a higher-order polynomial to compensate for the lack of the sliding window adaptivity in the end sections for real-time convolutional fitting. As a result, one can treat each fitting as the right-side end section and use the precalculated convolutional weights for efficient realtime smoothing and differentiation. Since the polynomial used in this case would have the correct frequency contents matching that of the experimental data being processed in a fixed window setting, the algorithm would be able to correctly detect the zerocrossing of the second derivative which corresponds to the inflection or end point of the titration curve (Figure 2D). However, every “outdated” point would also be “updated” with each slide of the filter window toward each newest point. This leads to additional flexibility as discussed before and results in excessive noise in the now “center” region.
CONCLUSION One cannot use the same order of polynomial fitting throughout the entire data set when using the nontruncating algorithm in a correct way. The signal frequency content in the end sections cannot be correctly matched using the same rigid polynomial appropriate for the center portion because of the lack of the sliding window adaptivity. Therefore, for the end sections one needs to increase the order of polynomial in order to increase the flexibility to match the real frequency contents of the signal being processed (so that one may obtain meaningful results). This is similar to splining techniques except convolutional least-squares fitting is used rather than direct point-to-point splining, which has no residual errors. Thus, one must use separate sets of polynomial fittings for the center and end sections: the order of polynomial used should be higher for the end sections than that used for the center section, in order to avoid effective end section truncations. Although in some cases the end sections obtained using the original nontruncating approach may be apparently adequate, they are in fact oversmoothed, as demonstrated in Figure 1. Two
modifications to the approach have been proposed and tested. For general postprocessing purposes, two different polynomial orders are needed: a lower-order one for the center section and a higherorder one for the end sections. For real-time feature detection purposes, such as end point detection in diffusional microtitration, a higher-order polynomial is required for effectively treating every newest m points as the right side end section. ACKNOWLEDGMENT This work was jointly funded by the NSF and the Whitaker Foundation (Cost Effective Health Care Technologies). It was also supported in part by the Elmer Lincoln Lindseth Chair of Biomedical Engineering at Case Western Reserve University. Received for review July 13, 1995. Accepted December 13, 1995.X AC950697K X
Abstract published in Advance ACS Abstracts, February 1, 1996.
Analytical Chemistry, Vol. 68, No. 6, March 15, 1996
1057
