Scipy filter coefficients and filtering in MATLAB

I generated filter coefficients with scipy's signal.butter to use in real time filtering in Arduino, and I'm validating the filter with a simple Matlab code that produces a dummy wave and produces a filtered wave from it.

I validated the scipy generated coefficients for a 1st order filter, but when doing 4th order the filter seems to become unstable, which makes no sense. I think it's a logic issue either on the math side, or the interpreting the coefficients side, but maybe it's a code issue.

Scipy generates 9 denominators (a0 to a8) and 5 numerators (I'm assuuming b0 to b4) but the matrix for b has the following format: [ 4.82121714e-13 0.00000000e+00 -1.92848686e-12 0.00000000e+00 2.89273028e-12 0.00000000e+00 -1.92848686e-12 0.00000000e+00 4.82121714e-13] In the first order validation I simply ignored the 0s, and assumed the third value to be b1, so I've been doing the same, assuming the non-0 values to be b0-b4. When running the code for 1st order the signal gets filtered (commented part of the code), but when doing it for 4th order it explodes (uncommented part of the code).

Matlab code

% --------- Data Señal de entrada
Fs = 20000;                 % frecuencia de muestreo (o cuantas muestras por segundo) (Hz)
T = 1/Fs;                   % periodo de muestreo (o delta t por muestreo) (s)
L = Fs;                     % Longitud de señal, L=Fs para ser 1 segundo
t = (0:L-1).*1/Fs;              % vector de tiempo (s)
f = 10;                     % frecuencia señal de entrada (Hz)
x = sin(2.*pi.*f.*t) + 0.1*sin(2*pi*1000*t) + 0.1*sin(2*pi*5000*t);     % señal de entrada con ruido

% -------- cutoff frecuencias, wc = 1/RC = 1/tao; wc = 2*pi*fc = 1/tao
fcLPF = 13;
fcHPF = f^2/fcLPF;
tao1 = 1/(2*pi*fcHPF);
tao2 = 1/(2*pi*fcLPF);

% -------- coeficientes, Band-pass 4to° orden
a0 = 1;

% a1 = -1.99832407;
% a2 = 0.99833393;

a1 = -7.99560326;
a2 = 27.96927189;
a3 = -55.90796275;
a4 = 69.84684949;
a5 = -55.84709416;
a6 = 27.90840316;
a7 = -7.96951656;
a8 = 0.99565219;

% b0 = 0.00083304;
% b1 = -0.00083304;

b0 = 4.82121714/10000000000000;
b1 = -1.92848686/1000000000000;
b2 = 2.89273028/1000000000000;
b3 = -1.92848686/1000000000000;
b4 = 4.82121714/10000000000000;



% ------- Algoritmo filtro Band-pass
y = zeros(1,length(t));

for n=9:length(t)
y(n) = - a1/a0*y(n-1) - a2/a0*y(n-2) - a3/a0*y(n-3) - a4/a0*y(n-4) - a5/a0*y(n-5) - a6/a0*y(n-6) - a7/a0*y(n-7) - a8/a0*y(n-8) + b0/a0*x(n) + b1/a0*x(n-1) + b2/a0*x(n-2) + b3/a0*x(n-3) + b4/a0*x(n-4);

% for n=3:length(t)
% y(n) = - a1/a0*y(n-1) - a2/a0*y(n-2) + b0/a0*x(n) + b1/a0*x(n-1);

end

% ------- figuras
figure;
hold on
plot(t,x,"k")
plot(t,y,"r","LineWidth",2)
ylabel("Amplitud de Señales")
xlabel("Tiempo [seg]")
legend("Señal Original","Señal Filtrada")
box on
hold off

I don't think it's an issue, but here's the python code I used, I'm running it on collab

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy import signal

bt,at= sp.signal.butter(4,[7.692307692307693,13],'bandpass',fs=20000)
b1,a1= sp.signal.butter(1,[7.692307692307693,13],'bandpass',fs=20000)

@Irreducible is correct. You have created filter coefficients that expose a numerical stability issue in the case of a 4th order Butterworth filter. See here for an explanation: https://en.wikipedia.org/wiki/Digital_filter#Direct_form_II

Basically, DFII is a highly optimized implementation of a time domain convolution, but since the poles are calculated first, the gain can go extremely high, extremely fast and lead to numerical instability, especially in the case of high Q filters (which yours is an example of).

Also, not sure why you see only 5 values for your bs, but scipy.signal.butter always returns the same number of coefficients for b and a.

As suggested, express the filter as second order sections and your filter will be fine:

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy import signal
fs=20000
bt,at= signal.butter(4,[7.692307692307693,13],'bandpass',fs=fs)
b1,a1= signal.butter(1,[7.692307692307693,13],'bandpass',fs=fs)

t = np.linspace(0, 1, fs, endpoint=False)
freq = 10
sqr_wave = signal.square(2 * np.pi * 10 * t)
plt.plot(t, sqr_wave,label='orig')
sqr_wave_1st_order = signal.lfilter(b1, a1, sqr_wave)
plt.plot(t, sqr_wave_1st_order,label='1st order')
sqr_wave_4th_order = signal.lfilter(bt, at, sqr_wave)
plt.plot(t, sqr_wave_4th_order,label='4th order ba')
sos_4th_order = signal.butter(4,[7.692307692307693,13],'bandpass',fs=fs,output='sos')
sqr_wave_sos = signal.sosfilt(sos_4th_order, sqr_wave)
plt.plot(t, sqr_wave_sos,label='4th order sos')
plt.ylim(-2, 2)
plt.legend()
plt.show()
print(bt)
print(at)

Output:

[ 4.82121714e-13  0.00000000e+00 -1.92848686e-12  0.00000000e+00
  2.89273028e-12  0.00000000e+00 -1.92848686e-12  0.00000000e+00
  4.82121714e-13]
[  1.          -7.99560326  27.96927189 -55.90796275  69.84684949
 -55.84709416  27.90840316  -7.96951656   0.99565219]
[ 0.00083304  0.         -0.00083304]

To implement your SOS filter, you basically take your set of biquads (output of ouptut='sos') and stack the filters back to back. Julius Smith has a succinct explanation: https://ccrma.stanford.edu/~jos/fp/Series_Second_Order_Sections.html.