0
$\begingroup$

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 ?

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Marcus Müller Sep 28 '17 at 9:05
  • $\begingroup$ @MarcusMüller : I added the link in the post. The paper does not say which filter is used. $\endgroup$ – Pouteri Sep 28 '17 at 9:07
  • 1
    $\begingroup$ @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. $\endgroup$ – Olli Niemitalo Sep 28 '17 at 12:57
  • 1
    $\begingroup$ @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 :) $\endgroup$ – Marcus Müller Sep 28 '17 at 14:21
  • $\begingroup$ @MarcusMüller I'm with Olli on this one. The context seems to imply taking a derivative, aka DC killer. $\endgroup$ – Cedron Dawg Jul 26 at 0:17
4
$\begingroup$

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)$.

Sinc and its derivative
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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – denis Nov 10 '17 at 15:39
  • $\begingroup$ @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. $\endgroup$ – Olli Niemitalo Nov 10 '17 at 23:18
0
$\begingroup$

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.
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.