Derivative filter in Python
Ideal derivative filter

Let $f(x)$ be a signal bandlimited to frequencies $(-\pi,\, \pi)$. Given $f(x)$ as input, the same $f(x)$ is given as output by a system that has as its impulse response the sinc function:

$$\operatorname{sinc}(x) = \left\{\begin{array}{ll}1&amp;\text{if }x = 0,\\
\frac{\sin(\pi x)}{\pi x}&amp;\text{otherwise.}\end{array}\right.\tag{1}$$

Taking the derivative $f'(x)$ of signal $f(x)$ is a linear time-invariant operation. By associativity, $f(x)$ is differentiated by a system that has as its impulse response the derivative of $\operatorname{sinc}(x):$

$$\operatorname{sinc}'(x) = \left\{\begin{array}{ll}0&amp;\text{if }x = 0,\\
\frac{\cos(\pi x)}{x} - \frac{\sin(\pi x)}{\pi x^2}&amp;\text{otherwise.}\end{array}\right.\tag{2}$$

Both of $\operatorname{sinc}(x)$ and $\operatorname{sinc}'(x)$ are bandlimited  to $(-\pi,\, \pi)$ and can thus be sampled at integer $x$ without aliasing. Sampling $\operatorname{sinc}'(x)$ at integer $x$ gives the ideal discrete-time impulse response for you to use in filtering samples of $f(x)$ at integer $x$ to obtain samples of $f'(x)$.

Figure 1. $\operatorname{sinc}(x)$ (red) and its derivative $\operatorname{sinc}'(x)$

Windowed filter

The infinitely long ideal impulse response can be multiplied by a window function to obtain a realizable impulse response. SciPy can calculate several types of window functions and do finite-impulse-response (FIR) filtering with an arbitrary impulse response by scipy.signal.convolve. You need to delay the impulse response to make it that of a causal system.

Least squares derivative filter

A least squares filter impulse response is obtained by sampling $\operatorname{sinc}'(x)$ from Eq. 2 symmetrically at integer $-N\le x\le N$ with $N$ determining the filter order. This includes the samples of the impulse response that have the largest absolute value and consequently have the largest contribution to the mean square error.

Least squares derivative filter for pre-oversampled data

For sampled data that is oversampled by a factor $\beta \ge 1$ compared to the sufficient sampling frequency of $1$, a least squares derivative filter is obtained by minimizing mean square error (MSE) or deviation of the Fourier transform of the impulse response from the Fourier transform $i \omega$ of derivation, over the bandwidth $\omega\in(-\pi,\,\pi)$ of sinc:

$$\begin{align}\mathrm{MSE} &amp;= \frac{1}{2\pi}∫_{-\pi}^{\pi}\left|iω - \sum_{n=1}^N \left(c_n e^{-niω/β} - c_n e^{niω/β}\right)\right|^2\\
&amp;= \frac{1}{\pi}∫_{0}^{\pi}\left(ω + 2 \sum_{n=1}^N c_n \sin\left(\frac{nω}{β}\right)\right)^2,\end{align}\tag{3}$$

where the unorthodox notation $e^{-iω/β}$ represents a one-sample delay at the oversampled sampling frequency or a delay of $1/β$ samples at the critical sampling frequency of the sinc, and $c_n$ is the value of the antisymmetric impulse response at discrete time index $n=1\ldots N$ with a complementary value $-c_{n}$ at time $-n$. The reason for the unorthodox notation is that it gives MS values that are comparable between different $\beta$. If you require coefficients for calculation of the derivative as if the sampling frequency was $1$, then divide each coefficient (solved in the following) by $\beta$. MSE is minimized when partial derivatives of MSE with respect to all $c_n$ are zero. At $\beta = 1$ we get the least squares filter discussed earlier. The least squares solutions and the corresponding impulse responses can be calculated by this Python script:

from sympy import *
for N in [1, 2, 3, 4]:  # &lt;------ number of non-zero coefs / 2, too large is too slow to solve
    omega, beta = symbols('omega beta', real=True)
    c = [Symbol('c_'+str(i + 1), real=True) for i in range(N)]
    MSE = integrate((omega + 2*sum([c[n]*sin((n + 1)*omega/beta) for n in range(N)]))**2, (omega, 0, pi))/pi
    for use_beta in [1, 1.5, 2, 4, 8]:  # &lt;------- Oversampling factor
        use_MSE = MSE.subs(beta, use_beta)
        sol = solve([diff(use_MSE, c[n]) for n in range(N)], [c[n] for n in range(N)])
        for n in range(N):
            print(str(-sol[c[N - n - 1]].evalf()), end=',')
        print('0', end=',')
        for n in range(N-1):
            print(str(sol[c[n]].evalf()), end=',')
        print(str(sol[c[N-1]].evalf()), end='  ')
        print('# N='+str(N)+', beta='+str(use_beta)+', MSE='+str(use_MSE.subs(sol).evalf()))
    print()


SymPy's solver will hang on some inputs, but manages to find the solutions for $N = 1\ldots4$ and $\beta = 1, 1.5, 2, 4, 8$:

-1.00000000000000,0,1.00000000000000  # N=1, beta=1, MSE=1.28986813369645
-1.13548530933771,0,1.13548530933771  # N=1, beta=1.5, MSE=0.178081981619031
-1.27323954473516,0,1.27323954473516  # N=1, beta=2, MSE=0.0475902571416442
-2.12680294939990,0,2.12680294939990  # N=1, beta=4, MSE=0.00251926307592025
-4.06211519814366,0,4.06211519814366  # N=1, beta=8, MSE=0.000151083662191318

0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000  # N=2, beta=1, MSE=0.789868133696453
0.330596796032270,-1.24876549439064,0,1.24876549439064,-0.330596796032270  # N=2, beta=1.5, MSE=0.0130608565500205
0.288948415598137,-1.51850657748724,0,1.51850657748724,-0.288948415598137  # N=2, beta=2, MSE=0.000919718098315458
0.382389721267460,-2.75841263168952,0,2.75841263168952,-0.382389721267460  # N=2, beta=4, MSE=2.58832725686194e-6
0.689919219439685,-5.37907341717520,0,5.37907341717520,-0.689919219439685  # N=2, beta=8, MSE=9.33496011568591e-9

-0.333333333333333,0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000,0.333333333333333  # N=3, beta=1, MSE=0.567645911474231
-0.116742660811357,0.423760333529959,-1.31069003715342,0,1.31069003715342,-0.423760333529959,0.116742660811357  # N=3, beta=1.5, MSE=0.00108046987847814
-0.0790631626589370,0.433017934639585,-1.64079658337694,0,1.64079658337694,-0.433017934639585,0.0790631626589370  # N=3, beta=2, MSE=2.01122717118286e-5
-0.0825922199951234,0.659247420212084,-3.07101401580036,0,3.07101401580036,-0.659247420212084,0.0825922199951234  # N=3, beta=4, MSE=3.01168290348385e-9
-0.140650397261751,1.22875551919474,-6.03556856432330,0,6.03556856432330,-1.22875551919474,0.140650397261751  # N=3, beta=8, MSE=6.53134493440793e-13

0.250000000000000,-0.333333333333333,0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000,0.333333333333333,-0.250000000000000  # N=4, beta=1, MSE=0.442645911474231
0.0443010798143126,-0.173517915867411,0.482863594553193,-1.34856865228067,0,1.34856865228067,-0.482863594553193,0.173517915867411,-0.0443010798143126  # N=4, beta=1.5, MSE=9.55482684708157e-5
0.0232132887981686,-0.144357365428612,0.528060960838853,-1.71358998904610,0,1.71358998904610,-0.528060960838853,0.144357365428612,-0.0232132887981686  # N=4, beta=2, MSE=4.70708311492379e-7
0.0191194948759654,-0.179273078273186,0.859556991552367,-3.25775951118075,0,3.25775951118075,-0.859556991552367,0.179273078273186,-0.0191194948759654  # N=4, beta=4, MSE=3.75015864227226e-12
0.0307242274505004,-0.317445566053152,1.62921582434694,-6.42899176483137,0,6.42899176483137,-1.62921582434694,0.317445566053152,-0.0307242274505004  # N=4, beta=8, MSE=4.88983660976190e-17


The filter with N=2, beta=2, MSE=0.000919718098315458 has the same number of taps as Rick Lyons' filter −3/16, 31/32, 0, −31/32, 3/16 that is slightly suboptimal in least squares sense with MSE=0.0009390870.