当前位置：首页 > news >正文

别再只会调库了！用Python从零推导二阶巴特沃斯滤波器的差分方程（附NumPy实现）

news 2026/7/3 3:51:31

从零构建二阶巴特沃斯滤波器：数学推导与NumPy实战

在信号处理领域，滤波器就像一位精准的调音师，能够从嘈杂的环境中提取出我们需要的频率成分。当你使用scipy.signal.butter函数时，是否曾好奇过这个黑盒子内部究竟发生了什么？本文将带你深入滤波器设计的数学核心，从模拟域到数字域，一步步推导出二阶巴特沃斯低通滤波器的差分方程，并用纯NumPy实现整个过程。

1. 巴特沃斯滤波器的数学本质

巴特沃斯滤波器的魅力在于它在通带内具有最大平坦的幅度响应。这种特性源于其传递函数的特殊构造方式。对于一个归一化的二阶巴特沃斯低通滤波器（截止频率为1 rad/s），其传递函数可以表示为：

H(s) = 1 / (s² + √2*s + 1)

这个看似简单的表达式背后隐藏着精妙的数学设计。让我们分解这个传递函数的组成部分：

分母多项式：决定了滤波器的极点位置
√2系数：确保通带最大平坦性的关键
单位常数项：使直流增益为1

在Python中，我们可以用符号计算来验证这个传递函数的频率响应：

import sympy as sp s = sp.symbols('s') H = 1/(s**2 + sp.sqrt(2)*s + 1)

当我们需要不同的截止频率时，需要进行去归一化处理。假设目标截止频率为ωc，只需将s替换为s/ωc：

H_denormalized = 1/((s/ωc)**2 + sp.sqrt(2)*(s/ωc) + 1)

2. 从模拟到数字：双线性变换的艺术

将模拟滤波器转换为数字滤波器有多种方法，而双线性变换因其保持稳定性等优点成为最常用的方法之一。这种变换的本质是将s域（拉普拉斯域）映射到z域（离散域）：

s = (2/T) * (1 - z⁻¹)/(1 + z⁻¹)

其中T是采样周期。这个变换不是简单的替换，它需要解决几个关键问题：

频率畸变问题：直接应用会导致频率响应扭曲
预畸变补偿：需要在设计阶段预先校正截止频率
代数运算复杂度：高阶滤波器会生成复杂的z域表达式

预畸变校正公式为：

ωa = (2/T) * tan(ωd*T/2)

其中ωa是模拟截止频率，ωd是数字截止频率。

在Python中实现预畸变：

import numpy as np def prewarp(digital_freq, sample_rate): T = 1.0 / sample_rate return (2/T) * np.tan(digital_freq * T / 2)

3. 差分方程的完整推导

现在我们将完整的推导过程分解为可操作的步骤：

从数字指标开始：确定数字截止频率fd和采样率fs
转换为角频率：ωd = 2πfd
预畸变处理：计算ωa = (2/T)tan(ωdT/2)
去归一化传递函数：将s替换为s/ωa
应用双线性变换：s = (2/T)(1-z⁻¹)/(1+z⁻¹)
整理为标准形式：得到关于z⁻¹的多项式比

经过一系列代数运算（建议使用符号计算工具验证），我们最终会得到标准形式的数字滤波器传递函数：

H(z) = (b0 + b1*z⁻¹ + b2*z⁻²) / (1 + a1*z⁻¹ + a2*z⁻²)

对应的差分方程为：

y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] - a2*y[n-2]

4. NumPy实现与验证

现在我们将上述数学推导转化为实际的Python代码。首先实现系数计算：

def design_butterworth_2order(fc, fs): """设计二阶巴特沃斯低通滤波器系数 参数: fc - 截止频率(Hz) fs - 采样率(Hz) 返回: b, a - 分子和分母系数 """ T = 1.0 / fs ωd = 2 * np.pi * fc # 数字截止频率(rad/s) ωa = (2/T) * np.tan(ωd * T/2) # 预畸变后的模拟频率 # 二阶巴特沃斯归一化系数 a0_norm = 1.0 a1_norm = np.sqrt(2) a2_norm = 1.0 # 去归一化 a0 = a0_norm * (ωa**2) a1 = a1_norm * ωa a2 = a2_norm # 双线性变换 b0 = a0 b1 = 2 * a0 b2 = a0 a0_prime = a0 + a1*(2/T) + a2*(4/(T**2)) a1_prime = 2*a0 - 2*a2*(4/(T**2)) a2_prime = a0 - a1*(2/T) + a2*(4/(T**2)) # 归一化系数 b = np.array([b0, b1, b2]) / a0_prime a = np.array([a0_prime, a1_prime, a2_prime]) / a0_prime return b, a[1:] # 返回b和a(去掉a0)

接下来实现滤波器应用函数：

def butterworth_filter(x, b, a): """应用IIR滤波器 参数: x - 输入信号 b - 分子系数 [b0, b1, b2] a - 分母系数 [a1, a2] 返回: y - 滤波后信号 """ y = np.zeros_like(x) # 初始化状态 x_1, x_2 = 0, 0 # x[n-1], x[n-2] y_1, y_2 = 0, 0 # y[n-1], y[n-2] for n in range(len(x)): y[n] = b[0]*x[n] + b[1]*x_1 + b[2]*x_2 - a[0]*y_1 - a[1]*y_2 # 更新状态 x_2, x_1 = x_1, x[n] y_2, y_1 = y_1, y[n] return y

为了验证我们的实现是否正确，我们可以与SciPy的结果进行对比：

from scipy import signal # 设计参数 fc = 10 # 截止频率10Hz fs = 100 # 采样率100Hz # 我们的实现 b, a = design_butterworth_2order(fc, fs) # SciPy实现 b_sp, a_sp = signal.butter(2, fc/(fs/2), btype='low') print("我们的系数:", b, a) print("SciPy系数:", b_sp, a_sp[1:]) # a_sp包含a0=1

如果实现正确，两者的系数应该非常接近（可能有微小差异源于计算顺序和精度处理）。

5. 优化实现与性能考量

直接实现差分方程虽然直观，但在实际应用中可能需要考虑更多因素：

数值稳定性：高阶滤波器可能需要采用级联结构
计算效率：利用NumPy的向量化操作可以大幅提升速度
状态保存：面向对象实现更适合实际应用

这里给出一个优化后的向量化实现：

class ButterworthFilter: def __init__(self, fc, fs): self.b, self.a = design_butterworth_2order(fc, fs) self.reset() def reset(self): self.x_buf = np.zeros(2) # 存储x[n-1], x[n-2] self.y_buf = np.zeros(2) # 存储y[n-1], y[n-2] def filter(self, x): y = np.zeros_like(x) for i in range(len(x)): y[i] = (self.b[0]*x[i] + self.b[1]*self.x_buf[0] + self.b[2]*self.x_buf[1] - self.a[0]*self.y_buf[0] - self.a[1]*self.y_buf[1]) # 更新缓冲区 self.x_buf[1], self.x_buf[0] = self.x_buf[0], x[i] self.y_buf[1], self.y_buf[0] = self.y_buf[0], y[i] return y def filter_vectorized(self, x): """向量化实现，效率更高""" y = np.zeros_like(x) # 使用lfilter的等效实现 y[0] = self.b[0] * x[0] y[1] = (self.b[0]*x[1] + self.b[1]*x[0] - self.a[0]*y[0]) for i in range(2, len(x)): y[i] = (self.b[0]*x[i] + self.b[1]*x[i-1] + self.b[2]*x[i-2] - self.a[0]*y[i-1] - self.a[1]*y[i-2]) return y

在实际项目中，还需要考虑以下工程实践问题：

初始瞬态问题：滤波器启动时需要处理初始条件
实时处理：对于流式数据需要维护状态
定点实现：在嵌入式系统中可能需要定点运算
多级串联：实现更陡峭的滚降特性

6. 应用实例：噪声信号滤波

让我们用一个完整的例子展示滤波器的实际效果。首先生成一个包含多个频率成分的测试信号：

# 生成测试信号 t = np.linspace(0, 1, fs) # 1秒时长 signal_freq = 2 # 信号频率2Hz noise_freq = 30 # 噪声频率30Hz clean = np.sin(2 * np.pi * signal_freq * t) noise = 0.5 * np.sin(2 * np.pi * noise_freq * t) x = clean + noise # 混合信号

然后应用我们设计的滤波器：

# 设计并应用滤波器 bw_filter = ButterworthFilter(fc=10, fs=fs) y = bw_filter.filter(x) # 绘制结果 import matplotlib.pyplot as plt plt.figure(figsize=(10, 6)) plt.plot(t, x, label='原始信号') plt.plot(t, clean, '--', label='干净信号') plt.plot(t, y, label='滤波后信号') plt.legend() plt.xlabel('时间(s)') plt.ylabel('幅值') plt.title('二阶巴特沃斯滤波器效果对比') plt.grid(True) plt.show()

通过这个例子，你可以直观地看到滤波器如何有效地抑制高频噪声，同时保留低频信号成分。调整截止频率fc，可以观察到滤波器特性的变化。