7
$\begingroup$

I read a paper about a brain-computer interface. In this paper the authors reported "each signal has been filtered with an 8-order band-pass Chebishev Type I filter which cut-off frequencies are 0.1 and 10 Hz and has been decimated according to the high cut-off frequency". I tried to design this filter with scipy:

import scipy.signal as signal
signal.cheby1(8,0.05,[0.1,10.0],btype='band',analog=0,output='ba')

The result was:

Warning: invalid value encountered in sqrt
(array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
    nan,  nan,  nan,  nan,  nan,  nan]), array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
    nan,  nan,  nan,  nan,  nan,  nan]))

I have no background in signal processing, so I actually don't know what I am doing. I don't know whether they used a IIR or FIR filter or whether I have to scale the cut-off frequencies or whether I'm using the wrong ripple. I hope you can help me.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
5
+50
$\begingroup$

The main issue with the example you gave is that the filter design function cheby1 is returning all NaNs, which isn't going to be a very good filter. The problem is how you're specifying the passband/stopband edge frequencies. This particular function is meant to emulate MATLAB's cheby1 function; the frequencies that you give it should be normalized, such that a value of 1 corresponds to half of the sample rate.

import scipy.signal as signal
fs = whatever_the_sample_rate_of_the_filter_input_is_going_to_be
signal.cheby1(8,0.05,[0.1/(fs/2),10.0/(fs/2)],btype='band',analog=0,output='ba')

I don't have SciPy handy, but that should at least correctly design the filter you want.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks, I will try that in the evening. My sampling rate is 240 Hz. Do you know how to apply the filter on the data? I found a formula on wikipedia (en.wikipedia.org/wiki/…), but this one gives me strange oscillating results (very big numbers...). So I guess I implemented something wrong. I implement this in C++, so I cannot call the function scipy.signal.lfilter. Btw. what will change when I change the ripple? $\endgroup$
    – alfa
    Apr 6, 2012 at 7:07
  • 2
    $\begingroup$ This is a numerically challenging filter, since you have poles very very close to the unit circle. You need to break the filter down into second order sections and apply those sequentially. $\endgroup$
    – Hilmar
    Apr 6, 2012 at 11:30
  • $\begingroup$ And how would I break the filter down into second order sections? I think Matlab has a function for this, but it is not available in scipy. $\endgroup$
    – alfa
    Apr 6, 2012 at 19:37
  • $\begingroup$ @alfa: no scipy does not have it, but I started translating octave code to python: gist.github.com/endolith/4525003 $\endgroup$
    – endolith
    Feb 13, 2013 at 4:27
2
$\begingroup$

Two cutoff frequencies typically means it's a bandpass filter with the highpass at 0.1Hz and the low pass at 10 Hz. The low pass cutoff (which is the higher of the two frequencies) determines by how much you can down sample. If your lowpass filter was infinitely steep, you could get away with a new sample rate of 20 Hz (twice the cutoff). Since it has limited steepness, you need to leave a guard band between the cutoff frequency and the new Nyquist frequency. How much you need depends on the order of the filter and how much aliasing noise you can tolerate.

In this specific example it seems they have down-sampled by a factor of 12 or thereabouts, which seems too aggressive to me and will likely result in a lot of aliasing noise.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for the explanations. I think the downsampling by a factor of 12 makes sense, because there is a lot of noise with high amplitude and frequency and the signal they were interested in had a very low frequency. I actually wanted to reproduce their results. Are there any practical tutorials (with code) on this kind of filters? $\endgroup$
    – alfa
    Mar 31, 2012 at 17:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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