0
$\begingroup$

I have a data_stream that I'm attempting to filter using two on-line butterworth bandpass filters to extract two frequency bands of interest. The original data_stream looks like this:

original

mean: -3.9
max: 1257.0
min: -504
std: 171.9

I then filter this signal for two different frequency bands of interest [0.04-0.149hz] and [0.15-0.4hz]:

# hz
sampling_rate = 130

# filters
nyq = sampling_rate / 2

# lf filter
lf_low = 0.04 # hz
lf_lowcut = lf_low / nyq
lf_high = 0.149999999 # hz
lf_highcut = lf_high / nyq
lf_b, lf_a = scipy.signal.butter(N = 5, Wn = [lf_lowcut, lf_highcut], btype = 'band', output = 'ba', fs = sampling_rate)
lf_z = scipy.signal.lfilter_zi(lf_b, lf_a)

# hf filter
hf_low = 0.15
hf_lowcut = hf_low / nyq
hf_high = 0.40
hf_highcut = hf_high / nyq
hf_b, hf_a = scipy.signal.butter(N = 5, Wn = [hf_lowcut, hf_highcut],  btype = 'band', output = 'ba', fs = sampling_rate)
hf_z = scipy.signal.lfilter_zi(hf_b, hf_a)

# lf filtered data
lf_filtered_data = []

# hf filtered data
hf_filtered_data = []

# filter
for i in range(len(data_stream)):
    
    # lf filter
    lf_out, lf_z = scipy.signal.lfilter(lf_b, lf_a, [data_stream[i]], zi = lf_z)
    lf_filtered_data.append(float(lf_out))

    # hf filter
    hf_out, hf_z = scipy.signal.lfilter(hf_b, hf_a, [data_stream[i]], zi = hf_z)
    hf_filtered_data.append(float(hf_out))

This results in filtered data that doesn't seem to be correct. The lower bandpass filter yields a signal with (1) many np.NaNs (2) filtered samples that increase exponentially in value from ~0 to 10^300 as we move along the signal.

>>>np.count_nonzero(np.isnan(np.asarray(lf_filtered_data)))
28413
>>>temp = np.asarray(lf_filtered_data)
>>>temp = temp[np.isfinite(np.asarray(lf_filtered_data))]
>>>temp
array([-1.20137942e-021, -1.59566811e-020, -7.63479728e-020, ...,
        7.72798739e+305,  8.03696256e+305,  8.35706676e+305])
>>>np.max(temp[0:500])
4.0634521811857574e-10
>>>np.max(temp[500:1000])
126.740622754278
>>>np.max(temp[1000:2000])
1.2445953926536413e+19
>>>np.max(temp[2000:3000])
2.2864052486603268e+36
...
>>>np.max(temp[-1000:])
8.357066763585627e+305

What am I doing wrong here? Thanks!

$\endgroup$

1 Answer 1

1
$\begingroup$

Filters with high order and very low cutoff frequencies are numerically tricky since the poles are extremely close to the unit circle. You don't have enough numerically precision to implement this in transfer function form.

You can probably do it in cascaded second order section instead, but you must design the filters in poles and zeros, not in numerator and denominator coefficients.

The choice of this filter is questionable. The settling time will be enormous: It will take over 2 and half minutes of signal to get an attenuation of 60 dB. Not sure what you are trying to do, but that doesn't look reasonable, even if you can implement it.

$\endgroup$
1
  • $\begingroup$ Thanks, appreciate the advice. I did end up trying to implement this using an online scipy sos filter, and it works much better. I'm fine with running the filter for a few minutes to get it settled. $\endgroup$ Apr 1, 2022 at 4:02

Your Answer

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

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