# Unexpected behaviour when using lfilter with initial filter delay multiple times

I want to filter a very long signal in smaller parts, therefor I am currently using scipys lfilter with an initial filter delay doing multiple iterations of:

(signal_part_filtered, filter_delay) = lfilter(b, a, signal_part, zi = filter_delay)


The first few iterarions are going as expected (original is blue; filtered is red):

But after a few iterations the inital filter delay keeps getting bigger:

I think the issue might be the limited precision of the float datatype, is there any way to solve this?

• What are your 'a' coefficients? – Matt L. Mar 12 '16 at 18:45
• I am using a butterworth bandpass filter (scipy's butter) with: a = [ 1. -9.37395837 39.56035793 -98.98354191 162.6115093 -183.27571567 143.52360875 -77.11103794 27.20283814 -5.68991959 0.53585937] – BStadlbauer Mar 13 '16 at 9:26
• Your filter is unstable, it has roots outside the unit circle. This means it can't work properly and will eventually result in overflow, regardless of precision. What are the filter specs? – Matt L. Mar 13 '16 at 9:45
• My samplerate is 160 samples/second; lowcut frequency is 0.1 Hz; highcut frequency is 5 Hz So I am calling b, a = butter(order, [low, high], btype='band') with: order = 5 low = 0.00125 high = 0.0625 – BStadlbauer Mar 13 '16 at 9:51

The problem is not in the filtering process but already at the design stage. Your specifications are very difficult to realize because your desired band is at very low frequencies. With these specs, the butter routine returns an unstable filter (i.e., with two pole pairs outside the unit circle of the complex plane).

One thing you could try is reduce the filter order. Always check the maximum pole radius by computing the (magnitude of the) roots of the denominator polynomial. It could be that the design turns out to be OK for a lower filter order. Of course, your filter will be less steep.

Another, better approach would be to downsample your signal before filtering, so your pass band is not anymore in such a low frequency range (compared to the sampling frequency). This will avoid numerical problems in the filter design routine.

Alternative suggestion: don't use polynomials [b,a] as filter representation but keep it in poles and zeroes [z,p,k]. Then execute the filter as cascaded second order section. This will insure the that the poles stay inside the unit circle. In Matlab that would look something like this

%% does not work
x0 = zeros(8192,1); x0(1) = 1; % digital delta function
[b,a] = butter(5,[0.00125 0.0625]);
y0 = filter(b,a,x0);
fprintf('Max of the impules response with polynomials = %g\n',max(abs(y0)));

%% does work
[z,p,k] = butter(5,[0.00125 0.0625]);
sos = zp2sos(z,p,k);
y1 = sosfilt(sos,x0);
fprintf('Max of the impules response with second order sections = %f\n',max(abs(y1)));