从零实现Chebyshev滤波器频率响应优化

从零构建Chebyshev滤波器：不只是“抄公式”的实战笔记

你有没有遇到过这种情况？系统里有个强干扰信号，紧贴着你的有用频带，Butterworth滤波器滚降太缓，压不住；换成高阶FIR吧，MCU直接卡死。这时候，工程师的工具箱里真正能打硬仗的，往往是那个名字拗口但性能彪悍的老将——Chebyshev Type I 滤波器。

它不追求通带“完美平坦”，而是聪明地牺牲一点小波动，换来过渡带近乎垂直的衰减。听起来像某种工程上的“以空间换时间”？没错，这正是它的精髓所在。

今天我们就抛开MATLAB的cheby1()函数，一行行手推数学、一步步写出代码，彻底搞懂如何在没有高级库支持的环境下，亲手打造一个稳定可用的Chebyshev滤波器。这不是理论复读机，而是一份可以贴到你嵌入式项目里的实战指南。

为什么是Chebyshev？当“平滑”不再是第一优先级

我们先来直面现实：在大多数教科书里，Butterworth被称为“最平坦”滤波器，听上去很美。但在真实世界中，“平坦”往往意味着代价——缓慢的滚降速度。

设想这样一个场景：你在设计一款便携式心电监测仪，采样率250Hz，想保留0.5–40Hz的心电信号，但环境中充斥着50Hz工频干扰。你用了一个4阶Butterworth低通，截止频率设为45Hz。结果呢？50Hz处的衰减可能只有15dB左右，噪声依然清晰可见。

这时候，如果你允许通带内有±0.5dB的轻微波动（人耳或生理信号处理几乎察觉不到），换成4阶Chebyshev Type I，同样的阶数下，50Hz处的衰减轻松突破30dB。这就是选择性带来的质变。

关键洞察：Chebyshev不是“更好”的Butterworth，它是另一种设计哲学——用可控的通带纹波换取极致的选择性。

核心机制拆解：切比雪夫多项式到底做了什么？

所有神奇都源于这个公式：

$$
|H(j\omega)|^2 = \frac{1}{1 + \varepsilon^2 C_n^2\left(\frac{\omega}{\omega_c}\right)}
$$

别被吓到，我们一帧一帧拆开看。

切比雪夫多项式：制造“可控震荡”的引擎

$ C_n(x) $ 是第n阶切比雪夫多项式，定义如下：

$$
C_n(x) =
\begin{cases}
\cos(n \cdot \arccos x), & |x| \leq 1 \
\cosh(n \cdot \text{arcosh } x), & |x| > 1
\end{cases}
$$

重点来了：
- 当 $ |\omega / \omega_c| \leq 1 $（也就是通带），$ x \in [-1,1] $，此时 $ C_n(x) $ 在 $[-1,1]$ 区间来回震荡（因为$\cos$函数振荡）；
- 而一旦进入阻带（$ |\omega / \omega_c| > 1 $），$ C_n(x) $ 变成双曲余弦，指数级增长 → 分母爆炸 → 增益急速下降。

所以你看，通带的“波纹”不是缺陷，而是设计使然；而阻带的陡峭衰减，正是来自双曲函数的“非线性放大”。

参数三剑客：你真正该调的是什么？

✅ 实战建议：
- 一般选 $ \delta_p = 0.1 \sim 1.0 $ dB，太大影响测量精度；
- 阶数优先用偶数（方便拆成二阶节SOS）；
- 不要盲目追求高阶，6~8阶已是多数MCU的极限。

手搓一个数字滤波器：从模拟极点到SOS系数

现在进入重头戏。我们要做的，是把上面的数学变成能在STM32或DSP上跑起来的差分方程。

第一步：生成模拟原型极点

Chebyshev滤波器的所有极点分布在s平面左半边的一个椭圆上。每个极点的位置由以下公式确定：

$$
\sigma_k = -\sinh\left(\frac{1}{n}\mathrm{arcsinh}\frac{1}{\varepsilon}\right)\cdot \sin\left(\frac{2k+1}{2n}\pi\right)
$$
$$
\omega_k = \cosh\left(\frac{1}{n}\mathrm{arcsinh}\frac{1}{\varepsilon}\right)\cdot \cos\left(\frac{2k+1}{2n}\pi\right)
$$

对应的复极点就是 $ p_k = \sigma_k + j\omega_k $。

这些极点最初是在归一化截止频率 $ \omega_c = 1 $ 下计算的，我们需要根据实际需求缩放。

第二步：双线性变换前的预畸变

数字滤波器不能直接照搬模拟频率，必须做预畸变防止“频率 warping”：

$$
\omega_c^{\text{prewarped}} = \tan\left( \pi \frac{f_c}{f_s} \right)
$$

然后将每个模拟极点乘上这个值：
p_analog = p_normalized * prewarp_wc

第三步：映射到z域

使用双线性变换关系 $ s = \frac{2f_s (z-1)}{z+1} $，反解出z域极点：

$$
z = \frac{1 + p/(f_s/2)}{1 - p/(f_s/2)}
$$

⚠️ 关键检查点：所有z域极点必须严格位于单位圆内！否则滤波器会发散。

第四步：组装成SOS结构

IIR滤波器直接用高阶传递函数极易因浮点误差失稳。正确做法是分解为多个二阶节（Biquad Sections）串联。

每两个共轭极点组成一个二阶节。若阶数为奇数，则单独留一个一阶节。

对于每个二阶节，其传递函数形式为：

$$
H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}
$$

分子初始可设为[1, 0, 0]，再通过归一化直流增益修正系数。

真·从零实现：Python版可移植核心代码

下面这段代码不含任何scipy.signal构造函数，完全手动推导，适合转写为C语言嵌入式版本。

import numpy as np def design_chebyshev_type_i(N, rp, fc, fs): """ 手动设计 Chebyshev Type I 低通滤波器 返回 SOS 系数矩阵 [nb_sections x 6] """ # 1. 预畸变截止频率 w_digital = 2 * fc / fs w_analog = np.tan(np.pi * w_digital / 2) # 2. 计算波纹因子 ε epsilon = np.sqrt(10**(0.1 * rp) - 1) # 3. 计算归一化模拟极点 poles_normalized = [] sinh_arg = np.arcsinh(1 / epsilon) / N for k in range(N): theta = np.pi * (2*k + 1) / (2*N) sigma = -np.sinh(sinh_arg) * np.sin(theta) omega = np.cosh(sinh_arg) * np.cos(theta) pole_norm = sigma + 1j * omega poles_normalized.append(pole_norm) # 缩放到目标截止频率 poles_analog = np.array(poles_normalized) * w_analog # 4. 双线性变换到z域 poles_z = (1 + poles_analog / (fs/2)) / (1 - poles_analog / (fs/2)) # 5. 分组为二阶节（SOS） sos = [] used = np.zeros(len(poles_z), dtype=bool) for i in range(len(poles_z)): if used[i]: continue # 寻找共轭对（或单实根） section_poles = [poles_z[i]] used[i] = True # 查找共轭 for j in range(i+1, len(poles_z)): if not used[j] and np.isclose(poles_z[i], np.conj(poles_z[j])): section_poles.append(poles_z[j]) used[j] = True break # 构造分母多项式 a(z) a = np.poly(section_poles) # 分子：默认无零点（全极点滤波器），但可加z=-1处零点改善高频响应 if len(section_poles) == 2: b = [1, 0, 0] # 可替换为 [1, -2, 1] 引入两个z=-1零点 else: b = [1, 0] # DC增益归一化 dc_b = np.polyval(b, 1) dc_a = np.polyval(a, 1) gain = abs(dc_b / dc_a) b = np.array(b) / gain # 补齐为6参数格式: [b0, b1, b2, a0, a1, a2] while len(b) < 3: b = np.append(b, 0) while len(a) < 3: a = np.append(a, 0) a[0] = 1.0 # 归一化a0 sos.append(np.concatenate([b, a])) return np.array(sos) # === 使用示例 === N = 4 rp = 0.5 # dB fc = 1000 # Hz fs = 8000 # Hz sos_custom = design_chebyshev_type_i(N, rp, fc, fs)

🔧 提示：此函数输出的sos结构与scipy.signal.sosfreqz兼容，可用于验证。

如何验证你没写错？对比才是硬道理

我们拿scipy的标准实现来做黄金参考：

from scipy import signal import matplotlib.pyplot as plt sos_scipy = signal.cheby1(N, rp, fc, btype='low', fs=fs, output='sos') # 计算频率响应 w, h_custom = signal.sosfreqz(sos_custom, worN=2000, fs=fs) _, h_scipy = signal.sosfreqz(sos_scipy, worN=2000, fs=fs) # 绘图 plt.figure(figsize=(10, 5)) plt.semilogx(w, 20*np.log10(np.abs(h_custom)+1e-10), label='Custom') plt.semilogx(w, 20*np.log10(np.abs(h_scipy)+1e-10), '--', label='SciPy') plt.grid(True, which='both') plt.xlabel('Frequency (Hz)') plt.ylabel('Magnitude (dB)') plt.legend() plt.title('Chebyshev Type I Frequency Response: Custom vs SciPy') plt.ylim(-60, 5) plt.show()

如果两条曲线基本重合，恭喜你，已经掌握了从数学到代码的完整闭环！

工程落地中的那些“坑”与应对策略

纸上得来终觉浅。当你真把这段代码烧进MCU时，以下几个问题大概率会出现：

❌ 坑点1：定点数导致极点越界，滤波器自激振荡

现象：输入静音，输出却持续振荡。
原因：系数量化后，某个极点从单位圆内跳到了外面。
对策：
- 使用Q31格式存储系数；
- 添加极点稳定性检查：确保|z| < 0.995；
- 或改用格型结构（Lattice Filter），天生更稳定。

❌ 坑点2：启动瞬间输出爆冲

现象：上电后第一秒输出极大值，甚至饱和。
原因：状态变量初值为0，但输入突然接入，相当于阶跃激励。
对策：
- 软启动：逐步增加增益（fade-in）；
- 或初始化状态为当前输入值的延迟副本。

❌ 坑点3：相位扭曲导致脉冲信号变形

现象：ECG的R波变得歪斜，音频听感“浑浊”。
原因：Chebyshev本质是非线性相位。
对策：
- 若允许离线处理，使用filtfilt类零相位滤波；
- 实时系统中可搭配全通网络做群延迟均衡。

实战应用场景举隅

场景1：EEG前端抗噪滤波

• 需求：提取0.5–30Hz脑电波，抑制50Hz及其谐波
• 方案：4阶Chebyshev I LPF，fc=35Hz,rp=0.5dB
• 效果：50Hz衰减 >40dB，CPU占用仅为同等性能FIR的1/10

场景2：SDR中的信道选择滤波

• 需求：在密集FM广播频段中分离单一电台
• 方案：6阶Chebyshev BP，中心频率100MHz±100kHz
• 优势：相比Butterworth节省2个阶数，降低FPGA资源消耗

场景3：工业传感器信号调理

• 需求：振动传感器输出含高频机械噪声
• 方案：定制化8阶Chebyshev，动态切换fc适应不同转速工况
• 特点：边缘设备自主重构，无需PC干预

写在最后：掌握原理，才能超越工具

今天我们走完了从切比雪夫多项式→模拟极点计算→双线性变换→SOS分解→代码实现与验证的完整链条。你会发现，一旦理解了背后的机制，哪怕没有MATLAB，也能在资源受限的嵌入式平台上精准控制每一个滤波行为。

更重要的是，这种“从零构建”的能力让你不再只是参数的使用者，而是成为系统性能的定义者。你可以：
- 动态调整波纹与阶数平衡实时性与抑制能力；
- 在FPGA中展开流水线提升吞吐；
- 结合机器学习在线优化滤波策略。

下次当你面对棘手的干扰问题时，希望你能想起这个“老派但高效”的武器——Chebyshev滤波器，以及你是如何亲手把它从数学符号变成运行在电路中的声音与数据的。

如果你正在做类似项目，欢迎在评论区分享你的实现细节或遇到的挑战，我们一起打磨这份“硬核”技能。