Scipy生成4阶巴特沃斯带通滤波系数在MATLAB验证不稳定问题

问题背景 #

使用scipy.signal.butter生成巴特沃斯带通滤波系数，用于Arduino实时滤波，通过MATLAB代码验证滤波效果时，1阶滤波表现正常，但4阶滤波出现结果发散的不稳定现象。

问题细节：

• Scipy生成的4阶带通滤波器b系数包含多个0值，手动提取非0值作为b0-b4使用
• a系数共9个（a0到a8），在MATLAB中实现4阶滤波时结果发散

重现代码 #

Python系数生成代码 #

import numpy as np
import scipy as sp
from scipy import signal

# 4阶带通滤波器，通带7.692~13Hz，采样率20000Hz
bt,at= sp.signal.butter(4,[7.692307692307693,13],'bandpass',fs=20000)
# 1阶带通滤波器，参数同上
b1,a1= sp.signal.butter(1,[7.692307692307693,13],'bandpass',fs=20000)

MATLAB滤波实现代码 #

% --------- 输入信号生成
Fs = 20000;                 % 采样率(Hz)
T = 1/Fs;                   % 采样周期(s)
L = Fs;                     % 信号长度，对应1秒
t = (0:L-1).*1/Fs;          % 时间向量(s)
f = 10;                     % 目标信号频率(Hz)
% 含噪声的输入信号：10Hz主信号 + 1000Hz、5000Hz噪声
x = sin(2.*pi.*f.*t) + 0.1*sin(2*pi*1000*t) + 0.1*sin(2*pi*5000*t);     

% -------- 截止频率计算（仅注释用，未实际参与系数生成）
fcLPF = 13;
fcHPF = f^2/fcLPF;
tao1 = 1/(2*pi*fcHPF);
tao2 = 1/(2*pi*fcLPF);

% -------- 4阶带通滤波器系数（手动提取并转换后的错误版本）
a0 = 1;
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 = 4.82121714/10000000000000;
b1 = -1.92848686/1000000000000;
b2 = 2.89273028/1000000000000;
b3 = -1.92848686/1000000000000;
b4 = 4.82121714/10000000000000;

% ------- 滤波算法实现
y = zeros(1,length(t));

% 4阶滤波循环（从第9个点开始）
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);
end

% ------- 绘图对比
figure;
hold on
plot(t,x,"k")
plot(t,y,"r","LineWidth",2)
ylabel("信号幅值")
xlabel("时间 [秒]")
legend("原始信号","滤波后信号")
box on
hold off

问题原因分析 #

1. 系数提取错误：scipy生成的n阶带通数字滤波器，分子(b)和分母(a)系数长度均为2n+1（4阶对应9个系数）。手动丢弃b系数中的0值，只保留5个非0值，导致滤波公式的分子项缺失，这是滤波发散的核心原因。
2. 数值精度损失：手动将科学计数法表示的系数转换为分数形式（如4.82121714e-13转为4.82121714/10000000000000），可能引入微小的精度误差，高阶滤波对这种误差的放大效应更明显，加剧不稳定。
3. 直接高阶IIR实现的稳定性问题：直接II型结构的高阶IIR滤波器容易因数值误差积累导致不稳定，尤其是在嵌入式系统（如Arduino）的有限精度环境下。

解决方法 #

1. 使用完整的滤波器系数 #

直接使用scipy输出的全部b和a系数，不要丢弃任何元素。例如4阶带通的b系数是9个值，在MATLAB实现时需要完整代入滤波公式：

% 正确的系数使用方式（直接复制scipy输出的bt和at）
b = [bt[0], bt[1], bt[2], bt[3], bt[4], bt[5], bt[6], bt[7], bt[8]];
a = [at[0], at[1], at[2], at[3], at[4], at[5], at[6], at[7], at[8]];

% 滤波循环修正
for n=9:length(t)
    y(n) = (-a(2)*y(n-1) -a(3)*y(n-2) -a(4)*y(n-3) -a(5)*y(n-4) ...
           -a(6)*y(n-5) -a(7)*y(n-6) -a(8)*y(n-7) -a(9)*y(n-8))/a(1) ...
           + (b(1)*x(n) +b(2)*x(n-1) +b(3)*x(n-2) +b(4)*x(n-3) +b(5)*x(n-4) ...
           +b(6)*x(n-5) +b(7)*x(n-6) +b(8)*x(n-7) +b(9)*x(n-8))/a(1);
end

2. 避免手动转换系数格式 #

直接复制scipy输出的原始浮点值（包括科学计数法形式）到MATLAB中，减少精度损失。

3. 采用二阶节级联实现（推荐用于Arduino） #

高阶巴特沃斯滤波器可以分解为多个二阶IIR节（biquad）级联，这种结构的数值稳定性远高于直接高阶实现，且更适合Arduino的资源限制。

使用scipy的signal.butter结合signal.zpk2sos生成二阶节系数：

# 生成二阶节系数
z, p, k = signal.butter(4, [7.692307692307693,13], 'bandpass', fs=20000, output='zpk')
sos = signal.zpk2sos(z, p, k)
print(sos)

每个二阶节的形式为：
$$y[n] = b_0x[n] + b_1x[n-1] + b_2x[n-2] - a_1y[n-1] - a_2y[n-2]$$
在Arduino中依次实现每个二阶节，前一个节的输出作为下一个节的输入即可。

4. 初始条件优化 #

直接II型实现时，初始的延迟单元（y(n-1)到y(n-8)）可以设为0，但如果仍有轻微不稳定，可以尝试用输入信号的前几个点初始化延迟单元，减少初始瞬态响应。

内容的提问来源于stack exchange，提问作者PinkyWho LNG Pinkywho4884