Scipy 滤波器系数和 MATLAB 中的滤波

Question

我使用 scipy 的 signal.butter 生成了滤波器系数，以便在 Arduino 中用于实时滤波，并且我使用一个简单的 Matlab 代码来验证滤波器，该代码生成一个虚拟波并从中生成一个滤波波。

我验证了 scipy 为一阶滤波器生成的系数，但是当执行四阶滤波器时，滤波器似乎变得不稳定，这没有意义。我认为这无论是在数学方面还是在解释系数方面都是一个逻辑问题，但也许这是一个代码问题。

Scipy 生成 9 个分母（a0 到 a8）和 5 个分子（我假设是 b0 到 b4），但 b 的矩阵具有以下格式： [ 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] 在第一顺序验证中，我只是忽略了 0，并假设第三个值是 b1，所以我一直在做同样的事情，假设非 0 值是 b0-b4。当运行第一阶代码时，信号被过滤（代码的注释部分），但是当执行第四阶时，信号会爆炸（代码的未注释部分）。

Matlab代码

% --------- 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

我不认为这是一个问题，但这是我使用的Python代码，我在协作上运行它

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)

Answer 1

@不可约是正确的。您创建的滤波器系数暴露了四阶巴特沃斯滤波器的数值稳定性问题。请参阅此处的解释：https://en.wikipedia.org/wiki/Digital_filter#Direct_form_II

基本上，DFII 是时域卷积的高度优化实现，但由于首先计算极点，增益可能会非常高、非常快，并导致数值不稳定，特别是在高 Q 值滤波器的情况下（您的滤波器是的一个例子）。

另外，不知道为什么你只看到

b

的 5 个值，但 scipy.signal.butter 总是为

b

和

a

返回相同数量的系数。

按照建议，将过滤器表示为二阶部分，您的过滤器就可以了：

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)

输出：

[ 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]

要实现 SOS 过滤器，您基本上需要使用双二阶集（ouptut='sos' 的输出）并将过滤器背靠背堆叠。朱利叶斯史密斯有一个简洁的解释：https://ccrma.stanford.edu/~jos/fp/Series_Second_Order_Sections.html.