# Derivative filter in Python

In Alaa Kharbouch, Ali Shoeb, John Guttag, Sydney S. Cash, An algorithm for seizure onset detection using intracranial EEG, Epilepsy & Behavior, Volume 22, Supplement 1, 2011 (section 2.1, 3rd paragraph) about EEG, the authors note that the spectral amplitude profile of a signal is inversely proportional to frequency. To correct for this trend, they propose to apply a derivative filter to the signal.

My question is: Is there a Python function which implements such a derivative filter ? Is the savgol_filter function from the Scipy module suited to this task ? If not, how could I design such a filter in Python ?

• Can you link to that paper? derivative only means that it's the derivative of some other filter, and as far as I can tell, this tells me nothing about what the frequency response of the filter would look like – Marcus Müller Sep 28 '17 at 9:05
• @MarcusMüller : I added the link in the post. The paper does not say which filter is used. – Pouteri Sep 28 '17 at 9:07
• @MarcusMüller Derivation has a magnitude frequency response $|\omega|$ and a phase frequency response $\frac{\pi}{2}$ for positive frequencies and $-\frac{\pi}{2}$ for negative frequencies. – Olli Niemitalo Sep 28 '17 at 12:57
• @OlliNiemitalo a filter that implements derivation, certainly, but a filter that is a derivative of another filter is something else, and because I've confused these two terms myself, I asked for original reference :) – Marcus Müller Sep 28 '17 at 14:21
• @MarcusMüller I'm with Olli on this one. The context seems to imply taking a derivative, aka DC killer. – Cedron Dawg Jul 26 '19 at 0:17

## 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&\text{if }x = 0,\\ \frac{\sin(\pi x)}{\pi x}&\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&\text{if }x = 0,\\ \frac{\cos(\pi x)}{x} - \frac{\sin(\pi x)}{\pi x^2}&\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} &= \frac{1}{2\pi}∫_{-\pi}^{\pi}\left|iω - \sum_{n=1}^N \left(c_n e^{-niω/β} - c_n e^{niω/β}\right)\right|^2\\ &= \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]:  # <------ 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]:  # <------- 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.

• Nice theory. But don't the very long tails make this difficult in practice ? If I'm not mistaken, the sum |sinc'(k)| at positive integers is 1 + 1/2 + 1/3 ... diverges. – denis Nov 10 '17 at 15:39
• @denis It diverges, but it only means that for some "pathological" bounded inputs like ..., -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, ... the output grows without bound as the window length is increased. It's the same kind of situation as with sinc interpolation, see question: Signal values we will 'miss' between sampling instances during sampling of band limited signals. – Olli Niemitalo Nov 10 '17 at 23:18

Here are some links and notes on various derivative filters, probably TL;DR .

For uniformly spaced data, use np.gradient ; this calculates central differences $(x_{i+1} - x_{i-1}) \, / \, h$, and fairs the edges nicely. For better filters, combine np.gradient s at different spacings, e.g.
$^4/_3 \ (x_{i+1} - x_{i-1}) \ - \ ^1/_3 \ (x_{i+2} - x_{i-2})$; see Finite_difference_coefficient.

Although it's common to use np.diff without thinking, one-sided differences $(x_{i+1} - x_i) \, / \, h$ amplify high-frequency noise — try [1 -1 1 -1 ...] . Use np.gradient instead.

You could use Savitzky-Golay. But you usually want to add constraints at various frequencies, e.g. derivfilter( [1 -1 1 -1 ...] ) = 0: playing around vs. off-the-shelf.

For very noisy data, there seem to be two main approaches:

1. filter first, e.g. with an order-2 Butterworth low-pass filter, before anything else
2. fit splines, which have smooth derivatives; see spline-fitted-1d-data on SO.