0
$\begingroup$

In a research paper (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
  • $\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
3
$\begingroup$

Let $f(x)$ be a signal band-limited 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.$$

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

Both of $\operatorname{sinc}(x)$ and $\operatorname{sinc}'(x)$ are band-limited 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)$

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.

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