Another page  on this site describes how to design "Chebyshev-optimal" filters to calculate numerical estimates of the first derivative from samples of a function. This page extends that approach to estimate the second derivative. These filters are also optimal in the sense of minimizing the maximum deviation from the "ideal design" specification. Despite promising results, the designs are not the best choice for all applications.
While the approach used for designing first-derivative approximators applies nicely, there are several details to take into account when applying the method to second-derivative approximators.
Filter response characteristic. Since the second derivative estimator is the equivalent of applying a first-derivative estimator twice, the frequency-domain response of the second-derivative filter will be the equivalent of two multiplications by
w and two phase shifts by
Π/2. Thus, in the critical low-frequency region, the response curve is an upward-curving quadratic curve. The combination of the two phase shifts is the equivalent of multiplying by a scalar factor
-1.0. This is a somewhat harder problem numerically because the upward curvature must be countered by downward curvature to return the response curve somewhat close to zero at high frequencies, avoiding excessive noise response. The filters tend to need a few more terms than the first-derivative approximators to achieve similar accuracy.
Basis functions. A second derivative filter will have even symmetry. Cosine waves have the desired properties. So, instead of using a sine-series approximation as in the first-derivative case, use a cosine series approximation.
No natural zero constraint. The second derivative of a constant sequence is zero, much like the first derivative. However, cosine "basis functions" do not have a natural property of passing through zero. Thus, a constraint that the combination of cosine functions must equal 0 at 0 frequency must be included, with rigorous tolerance on this requirement. Since the value of a cosine function is 1.0 when evaluated at 0, this constraint is equivalent to saying that the sum of the filter coefficients must be equal to 0.
Basis size. The number of terms in each of the basis vectors of matrix
N/2 - Nt terms rather than
N/2 - 1 - Nt terms, to account for the added constraint at location 0.
After taking these adjustments into account, the solution method is fundamentally the same as for the first-derivative case. Full details of the method are provided on the page that covers the first-derivative case .
Here is the result of a typical design (see the 15-term design below).
An m-file design script for Octave shows the numerical details of constructing the matrix terms, calculating the solution, and plotting the results for the second-derivative filter design. If you wish to experiment with this script, download a copy of the minmax2.m file. The Octave solution is based on the GLPK package, but you can use the equivalent solver from the Optimization Toolbox for Matlab .
To see the practical benefits of the min-max designs, the 15-term filter (as tabulated below) is compared to the 9-term "Robust Low Noise" design from Pavel Holoborodko's site. While this comparison seems inherently unfair — the min-max design is allowed 60% more parameters to work with — this choice is pretty much forced.
The ninth-order "Robust Low Noise" designs is the highest order filter listed.
The "Robust Low Noise" designs deliver worse accuracy as the filter order increases, as you can clearly see from the frequency response plots for estimators of various orders — more parameters only makes the comparison more unfair.
First, here is the performance delivered by the RLN method in the presence of white random noise at 5% of the signal level, at a frequency of about 8% of the sampling rate (approximately 12 terms span a complete waveform cycle). A mathematically perfect result in green is plotted on top of the calculated results in blue.
Next, here is the performance delivered by the min-max method, same level of random noise, same signal. (The apparent different number of terms is only an edge effect artifact, nothing of significance. When a longer filter is convolved with the same length data block, there are more invalid values at the ends where the filter "overlaps" beyond the ends of the data block.)
At first glance, it might seem like the min-max design is undisputed winner — but unfortunately, the choice is not that clear. While the min-max design tracks changing curvature and higher frequencies better, in situations where the curvature is small and changing slowly, a "Robust Low Noise" design produces superior results. This is not yet fully understood, but it is consistent with the contrasting design strategies. The "RLN" approach concentrates maximum effort at matching low frequencies accurately. In contrast, the "min-max" design tolerates a small absolute error at every frequency, and this can result in a significant "relative error" when the size of the desired result is also supposed to be small. For now, the min-max designs are not a complete and general solution.
The following table provides some design results that you can use directly without resorting to the design script. For each target number of filter terms, the design parameters, the resulting performance specs, and filter coefficients are listed. If you are not sure about where to start, try the 15-term filter design with its
0.08·Fs guaranteed accurate band. If you need more economy or better performance, switch to one of the other designs.
|11-term — 5% accurate — 8% useful|
|Design: pass: 0.05 transition: 0.185 sensitivity: 100.0|
|pass ripple: 0.001 stopband edge: 0.24 stopband ripple: 0.10|
|13-term — 6% accurate — 9% useful|
|Design: pass: 0.08 transition: 0.135 sensitivity: 100.0|
|pass ripple: 0.001 stopband edge: 0.21 stop ripple: 0.104|
|15-term — 8% accurate — 11% useful|
|Design: pass: 0.08 transition: 0.14 sensitivity: 150|
|pass ripple: 0.00075 stopband edge: 0.22 stopband ripple: 0.1|
|17-term — 10% accurate — 13% useful|
|Design: pass: 0.10 transition: 0.16 sensitivity: 420|
|pass ripple: 0.0002 stopband edge: 0.26 stopband ripple: 0.082|
 The page on this site about first-derivative filter approximation has the details about to set up and solve the approximation problem.
 The "robust low noise estimators" are described at Pavel Holoborodko's web site, http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/.
 You can find more information about the Optimization Toolkit at
Site: RidgeRat's Barely Operating Site http://home.earthlink.net/~ltrammell Author: Larry Trammell Created: Dec 19, 2011 Revised: Jan 28, 2012 Status: Experimental Contact: NOSPAM ltramme1476 AT earthlink DOT net NOSPAM Restrictions: This information is public.