0
$\begingroup$

I'm designing a quite unusual (in terms of signal analysis) EEG experiment, which forces me to design my own filter. I have never done this before, so kindly please check my work and suggest corrections.

The analysis consists of cutting a narrow frequency band (2Hz) from 6-second EEG recordings(4*) (sampling frequency 256 Hz). I expect that after such filtering, I will receive a "somehow oscillated" signal (not in a strict mathematical sense, just with notable local maxima (peaks) and local minima (troughs))(5*). Now I'm interested just at the very end of this 6-sec window, the last dozen or so data points. I want to determine (by the simple algorithm provided by myself) if the signal ends "near the peak" or "near the trough" or in either of the above. So, all 6-sec windows are compressed to one of the three categories. You can think of this as a very rough estimate of the phase at the end of the window.

The tricky part is that I need to analyze overlapping frequency bands ( 4-6 Hz, 4.25 - 6.25 Hz, 4.5 - 6.5 Hz, etc.), so it will be great to have an extremely steep roll-off even at the price of amplitude disturbance (most of this information will be lost anyway) - that's why I came up with an elliptic IIR filter. Zero-phase shifting is also crucial so, forward-backward (filtfilt) filtering is applied. However, the devil is in the details (filter order, maximum ripple, minimum attenuation), and I must honestly admit that I determined these parameters by trial and error.

from scipy import signal
from mne.viz import plot_filter

sfreq = 256.
f_s = 6.
f_p = 8.
flim = (f_s - 1.0, f_p + 1.0)  # limits for plotting

nyq = sfreq / 2.  # the Nyquist frequency is half our sample rate
freq = [0, f_s, f_s, f_p, f_p, nyq]
gain = [0, 0, 1, 1, 0, 0]

ftype = "ellip"
order = 8
sos = signal.iirfilter(order, [f_s / nyq, f_p/nyq],
                       btype='bandpass', rp=5, rs=35, ftype=ftype, output='sos')
plot_filter(dict(sos=sos), sfreq, freq, gain,
            f'{ftype} order={order}', flim=flim, compensate=True)

Filters Visualisations:

Code results

So my questions are as follow:

  1. Is my filter "reasonable"?
  2. Can I do better?
  3. Can you recommend to me some learning materials for a filter design.
  4. (*) This is obviously cut from continuous EEG signal - can it be just cut (boxcar window), or should I apply something different (eg. Hamming window)?
  5. (*) Is my expectation justified?
$\endgroup$
0
$\begingroup$

Is my filter "reasonable"?

I would say no. The time domain ringing time of your filter far exceeds the length of your signal, there isn't enough signal for the output to stabilize.

Can I do better?

Difficult to say. Filter design involves complicated tradeoffs between many different factors and the "best" really depends on the specific needs & sensitivities of your application. You have partially explained what you want to do but not in enough details to give a recommendation. I'm skeptical about "looking at the phase" approach. Signal phase has only meaning if you have a well defined time reference for $t=0$, which most signals don't.

Can you recommend to me some learning materials for a filter design.

I personally like that one https://ccrma.stanford.edu/~jos/filters/ but I'm also mostly and audio guy.

(*) This is obviously cut from continuous EEG signal - can it be just cut (boxcar window), or should I apply something different (eg. Hamming window)?

Depends a bit on your application details. Ideally your filter "ringing" is significantly shorter than your signal. Then you can just take a "large" chunk of signal, filter it and cut off the edges to remove the filter transients.

(*) Is my expectation justified?

Can't say without looking at your signal first. A sufficiently narrow band signal will look like an amplitude and/or phase modulated sine wave. However, if your original signal is very noise-like, you'd have to make it VERY narrow to look like that.

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