Questionable results of on-line butterworth bandpass filtering using scipy.signal.lfilter

Question

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:

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!

Answer 1

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.