# Frequency response of numerical differentiation by polynomial interpolation / finite difference

One can use polynomial interpolation (or finite difference) to do numerical differentiation. However, there seems to be a surprising lack of interest in obtaining frequency response of this numerical differentiation method. Is there any reference I can look up to, or can anyone provide what frequency response would be? I want to obtain frequency response of any filter order (any number of filter taps).

By higher order, I mean higher-order differences as in https://en.wikipedia.org/wiki/Finite_difference#Higher-order_differences or in case of polynomial interpolation, polynomial degree.

• Could you define mathematically which type of higher order numerical differentiation you refer to? – Matt L. Dec 31 '18 at 9:19
• @MattL. Higher-order differences as in en.wikipedia.org/wiki/… or in case of polynomial interpolation, polynomial degree. – user8059 Dec 31 '18 at 9:24
• The freqz() function in Matlab/Octave can do this for you numerically. – Andy Walls Dec 31 '18 at 12:40

Polynomials are not integrable, hence their Fourier interpretation is complicated. However, locally, there are a few works. Under least-squares, polynomial interpolation's linearity works well with frequency analysis.

The class of Savitsky-Golay filters performs least-squares polynomial interpolation, and yields different derivative orders. As they seemed to be little known in signal processing, R. W. Schafer dedicated them What Is a Savitzky-Golay Filter? (public version) in the Lecture Notes section of IEEE Signal Processing Magazine in 2011. There, he especially covers their properties, including frequency responses. Analysis was also performed in Smoothing by Savitzky-Golay and Legendre filters by P. O. Persson and G. Strang in 2003. One can combine derivative accuracy and noise shaping or frequency response fitting.

More recently, there were derivations with IIR filters and fractional differentiation, for instance Maximally Flat IIR Smoothers With Repeated Poles and a Prescribed Delay, and even non-linear versions.

In image processing, there has been a lot of works as well, I can dig into my reference database if needed.

For finite differences, as suggested by Andy, you can use MATLAB/Octave's freqz, but with some tricks to get the phase frequency response right. I'm using Wikipedia's definitions of higher-order finite differences. In Octave:

graphics_toolkit("gnuplot");
n = 2; #finite differences derivative order
s = 0; #-1 for forward, 0 for central, 1 for backward
b = [1];
a = [1];
for k=1:n,
b = conv(b, [1, 0, -1]);
a = conv(a, horzcat(zeros(1, 1-s), [1]));
end;
w = [(0:512)/512*pi-pi/2];
h = freqz(b, a, w);
freqz_plot(w*2/pi, h, true);


Feel free to change n and s. b will contain the filter coefficients, diluted by zeros to trick freqz to calculate frequency responses when taps may be at half-integer times, and a will have a 1 to mark the tap at time zero in the impulse response. For the above:

b =
1   0  -2   0   1

a =
0   0   1


These are coefficients of the numerator and denominator polynomials of $$z^{-1}$$ of the transfer function, respectively. For forward and central finite differences the filter is acausal. The difference to the causal backward finite differences filter is that the impulse response has been shifted backwards $$m$$ samples by multiplying the transfer function by $$z^m$$ with positive $$m.$$ This introduces negative powers of $$z^{-1}$$ such as $$(z^{-1})^m$$ in the numerator. This can be made valid input to freqz by multiplying both the numerator and denominator by $$z^{-m}$$ to ensure both are polynomials of $$z^{-1}$$. Effectively, the script above does that.

Figure 1. freqz_plot output for central 2nd derivative finite differences.

The magnitude frequency response will be equal for all: forward, backward, central.