1. 项目缘起：为什么我们需要三参数陷波滤波器？

在信号处理、音频工程、振动控制以及通信系统设计中，我们常常会遇到一个棘手的问题：如何精准地、可调节地滤除一个特定频率的干扰信号，同时尽可能少地影响该频率附近的其他有用信号？这就是陷波滤波器（Notch Filter）的用武之地。标准的二阶陷波滤波器，其传递函数通常由中心频率和品质因数（Q值）两个参数决定。它像一个精准的手术刀，能在频谱上“挖”掉一个很窄的坑。但在实际工程中，尤其是自适应滤波、主动噪声控制或需要在线调整滤波器特性的场景里，我们常常希望这个“坑”的深度也是可控的。换句话说，我们不仅想决定“坑”的位置（频率）和宽度（Q值），还想决定“坑”到底有多深。这就是三参数陷波滤波器诞生的背景。

我最初接触这个需求是在一个电机驱动器的谐波抑制项目中。电机运行时会产生特定频率的电磁噪声，这个噪声频率会随着转速变化。我们需要一个滤波器能实时跟踪并抑制这个变化的噪声。标准陷波器虽然能跟踪频率，但其在中心频率处的衰减是固定的（理论上无穷大，实际受限于数值精度），这有时会“矫枉过正”，把一些有用的谐波成分也过度削弱了，影响了系统的动态响应。这时，一个具有独立衰减深度控制参数的三参数陷波器就显得非常必要。它允许我们在抑制干扰和保留信号完整性之间取得一个更灵活的平衡。

网络上关于标准陷波器的资料很多，但系统阐述三参数陷波器、特别是其从连续域（s域）到离散域（z域）完整推导过程的中文资料相对零散。很多MATLAB实现也只是直接调用iirnotch函数，对其背后的第三个参数“深度”如何融入传递函数，以及离散化时需要注意什么，往往语焉不详。本文将从一个一线工程师的视角，手把手推导三参数陷波滤波器的离散化过程，并给出可直接复现的MATLAB代码，重点解释第三个参数的意义和离散化带来的影响。

2. 三参数陷波滤波器的核心：传递函数剖析

要理解离散化，必须先吃透其在连续时间域（s域）的模型。一个三参数陷波滤波器的标准传递函数形式如下：

[ H(s) = \frac{s^2 + \omega_0^2}{s^2 + \frac{\omega_0}{Q}s + \omega_0^2} ]

这是经典的双二阶陷波形式。其中：

• (\omega_0 = 2\pi f_0)：是陷波的中心角频率（rad/s），对应需要滤除的干扰频率 (f_0)。
• (Q)：品质因数，决定了陷波的宽度。Q值越高，陷波越窄；Q值越低，陷波越宽。

这个传递函数在 (s = j\omega_0) 时，分子为零，分母不为零（假设Q有限），因此理论上在 (\omega_0) 处的增益为零，即无限衰减。但它的衰减深度是固定的、最大的。为了引入深度控制，我们需要修改传递函数，使其在中心频率处的增益不为零，而是一个我们可以设定的值 (\epsilon)（其中 (0 \le \epsilon < 1)）。(\epsilon=0) 代表完全抑制，(\epsilon=1) 则相当于全通（无任何滤波效果）。

一种常见且数学上优雅的引入方式是将分子项进行加权：

[ H_{3p}(s) = \frac{s^2 + \beta \cdot \frac{\omega_0}{Q}s + \omega_0^2}{s^2 + \frac{\omega_0}{Q}s + \omega_0^2} ]

这里，我们引入了一个新的参数 (\beta)。注意，分子中多了一项 (\beta \cdot \frac{\omega_0}{Q}s)。这个项就是控制深度的关键。

为什么是 (\beta)？让我们来分析一下：

1. 当 (\beta = 1) 时，分子和分母完全相同，(H_{3p}(s) = 1)，滤波器变成一个全通网络，没有任何滤波效果。这对应深度最浅（无衰减）的情况。
2. 当 (\beta = 0) 时，分子变回 (s^2 + \omega_0^2)，这就是标准的二阶陷波滤波器，在中心频率处衰减最大。这对应深度最深的情况。
3. 当 (0 < \beta < 1) 时，分子在 (s = j\omega_0) 处的值不再为零。我们可以计算出滤波器在中心频率 (\omega_0) 处的幅值响应 (|H_{3p}(j\omega_0)|)： [ H_{3p}(j\omega_0) = \frac{(j\omega_0)^2 + \beta \cdot \frac{\omega_0}{Q}(j\omega_0) + \omega_0^2}{(j\omega_0)^2 + \frac{\omega_0}{Q}(j\omega_0) + \omega_0^2} = \frac{-\omega_0^2 + j\beta \frac{\omega_0^2}{Q} + \omega_0^2}{-\omega_0^2 + j\frac{\omega_0^2}{Q} + \omega_0^2} = \frac{j\beta \frac{\omega_0^2}{Q}}{j\frac{\omega_0^2}{Q}} = \beta ] 看，结果非常简洁！在中心频率 (\omega_0) 处，滤波器的增益正好等于参数 (\beta)。因此，我们可以通过直接设定 (\beta) 的值，来精确控制滤波器在陷波中心处的衰减深度。例如，设定 (\beta = 0.1)，意味着在 (f_0) 处的信号幅度会被衰减到原来的10%（即-20 dB）。这个关系清晰直观，是这种形式的三参数陷波器被广泛采用的主要原因。

注意：参数 (\beta) 与之前提到的 (\epsilon) 是同一个概念，即 (\epsilon = \beta)。有些文献也可能用 (g) 或depth表示。在工程实现中，我们更关心其物理意义：中心频率处的增益。

3. 从连续到离散：双线性变换法详解

我们设计的 (H_{3p}(s)) 是模拟（连续时间）滤波器。要在数字系统（如DSP、FPGA或MATLAB/Simulink）中实现它，必须将其离散化，转化为适用于数字处理的差分方程形式，即z域传递函数。双线性变换（Bilinear Transform）是IIR滤波器离散化最常用、最稳定的方法之一，尤其适用于频率响应匹配。

双线性变换的核心公式是： [ s = \frac{2}{T} \cdot \frac{1 - z^{-1}}{1 + z^{-1}} ] 其中，(T) 是数字系统的采样周期，(f_s = 1/T) 是采样频率。

这个公式的妙处在于，它将s平面的左半平面（稳定区域）唯一地映射到z平面的单位圆内（稳定区域），保证了稳定性。但代价是引入了频率扭曲（Frequency Warping）：s域的模拟频率 (\omega_a) 和z域的数字频率 (\omega_d) 之间是非线性关系 (\omega_a = \frac{2}{T} \tan(\frac{\omega_d T}{2}))。

这对我们设计陷波器意味着什么？我们不能简单地把设计好的模拟中心频率 (f_0) 直接套用。如果我们希望离散后的数字滤波器在数字频率 (\omega_d = 2\pi f_0 / f_s) 处出现陷波，那么我们在设计模拟原型滤波器时，必须使用一个“预畸变”（Pre-warped）的模拟频率 (\omega_0')： [ \omega_0' = \frac{2}{T} \tan\left(\frac{\omega_d T}{2}\right) = \frac{2}{T} \tan\left(\frac{\pi f_0}{f_s}\right) ]这是整个离散化过程中最关键、也最容易出错的一步。很多初学者实现的滤波器频率不准，问题就出在忽略了预畸变。

我们的离散化步骤如下：

1. 确定数字指标：给定采样频率 (f_s)，期望陷波的数字中心频率 (f_d = f_0)，品质因数 (Q)，陷波深度 (\beta)。
2. 预畸变计算：计算预畸变后的模拟角频率 (\omega_0' = 2 \pi f_0' = \frac{2}{T} \tan(\pi f_0 / f_s))。
3. 构造模拟传递函数：将 (\omega_0') 代入三参数模拟传递函数： [ H_{3p}(s) = \frac{s^2 + \beta \cdot \frac{\omega_0'}{Q}s + (\omega_0')^2}{s^2 + \frac{\omega_0'}{Q}s + (\omega_0')^2} ]
4. 应用双线性变换：将 (s = \frac{2}{T} \cdot \frac{1 - z^{-1}}{1 + z^{-1}}) 代入上式。
5. 整理为z域标准形式：经过繁琐但必要的代数运算，将结果整理成关于 (z^{-1}) 的有理多项式形式： [ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{a_0 + a_1 z^{-1} + a_2 z^{-2}} ] 通常会将分母首项系数 (a_0) 归一化为1，得到： [ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-1}}{1 + a_1 z^{-1} + a_2 z^{-2}} ] 这里的系数 (b_0, b_1, b_2, a_1, a_2) 就是我们实现差分方程所需的系数。

4. 系数推导：一步步的手算过程

让我们把第3步的代数运算展开。为了简化书写，令： [ K = \frac{2}{T}, \quad \Omega_0 = \omega_0' = K \tan\left(\frac{\pi f_0}{f_s}\right), \quad \alpha = \frac{\Omega_0}{Q} ] 则模拟传递函数为： [ H(s) = \frac{s^2 + \beta \alpha s + \Omega_0^2}{s^2 + \alpha s + \Omega_0^2} ]

代入双线性变换 (s = K \frac{1 - z^{-1}}{1 + z^{-1}})：

分子 (N_s)： [ \begin{aligned} N_s &= \left[K \frac{1 - z^{-1}}{1 + z^{-1}}\right]^2 + \beta \alpha \left[K \frac{1 - z^{-1}}{1 + z^{-1}}\right] + \Omega_0^2 \ &= \frac{K^2(1 - z^{-1})^2}{(1 + z^{-1})^2} + \frac{\beta \alpha K (1 - z^{-1})}{1 + z^{-1}} + \Omega_0^2 \ &= \frac{K^2(1 - 2z^{-1} + z^{-2}) + \beta \alpha K (1 - z^{-1})(1 + z^{-1}) + \Omega_0^2 (1 + z^{-1})^2}{(1 + z^{-1})^2} \ &= \frac{K^2(1 - 2z^{-1} + z^{-2}) + \beta \alpha K (1 - z^{-2}) + \Omega_0^2 (1 + 2z^{-1} + z^{-2})}{(1 + z^{-1})^2} \end{aligned} ]

分母 (D_s)： [ \begin{aligned} D_s &= \left[K \frac{1 - z^{-1}}{1 + z^{-1}}\right]^2 + \alpha \left[K \frac{1 - z^{-1}}{1 + z^{-1}}\right] + \Omega_0^2 \ &= \frac{K^2(1 - 2z^{-1} + z^{-2}) + \alpha K (1 - z^{-1})(1 + z^{-1}) + \Omega_0^2 (1 + z^{-1})^2}{(1 + z^{-1})^2} \ &= \frac{K^2(1 - 2z^{-1} + z^{-2}) + \alpha K (1 - z^{-2}) + \Omega_0^2 (1 + 2z^{-1} + z^{-2})}{(1 + z^{-1})^2} \end{aligned} ]

因此，z域传递函数为： [ H(z) = \frac{N_s}{D_s} = \frac{K^2(1 - 2z^{-1} + z^{-2}) + \beta \alpha K (1 - z^{-2}) + \Omega_0^2 (1 + 2z^{-1} + z^{-2})}{K^2(1 - 2z^{-1} + z^{-2}) + \alpha K (1 - z^{-2}) + \Omega_0^2 (1 + 2z^{-1} + z^{-2})} ]

现在，合并分子分母中 (z^{-0}, z^{-1}, z^{-2}) 的系数。这是一个纯体力活，但必须仔细。

令公共分母为 (A = K^2 + \alpha K + \Omega_0^2)。实际上，我们通过整理分子分母的同类项，可以得到一组对称的系数。经过整理（过程略，建议读者自行推导一遍以加深理解），我们得到归一化后（即令分母常数项为1）的差分方程系数：

定义： [ D = K^2 + \alpha K + \Omega_0^2 ]

则分子系数（(b_i)）： [ \begin{aligned} b_0 &= (K^2 + \beta \alpha K + \Omega_0^2) / D \ b_1 &= (-2K^2 + 2\Omega_0^2) / D \ b_2 &= (K^2 - \beta \alpha K + \Omega_0^2) / D \end{aligned} ]

分母系数（(a_i)， 注意 (a_0=1)）： [ \begin{aligned} a_1 &= (-2K^2 + 2\Omega_0^2) / D \ a_2 &= (K^2 - \alpha K + \Omega_0^2) / D \end{aligned} ]

观察一下，你会发现 (b_1 = a_1)。这是这种特定形式的三参数陷波器经过双线性变换后的一个性质。同时，当 (\beta = 1) 时，(b_0=1, b_1=a_1, b_2=a_2)，滤波器确为全通；当 (\beta=0) 时，(b_0 = (K^2 + \Omega_0^2)/D), (b_2 = (K^2 + \Omega_0^2)/D)，变为标准陷波。

对应的差分方程为： [ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2] ] 其中 (x[n]) 是输入信号，(y[n]) 是输出信号。

5. MATLAB实现：从理论到代码

理论推导完成，接下来就是用MATLAB将其实现。我们将编写一个函数three_param_notch，输入设计参数，返回滤波器系数，并绘制频率响应图进行验证。

function [b, a] = three_param_notch(f0, Q, beta, fs) % 三参数陷波滤波器设计 (基于双线性变换) % 输入: % f0 - 陷波中心频率 (Hz) % Q - 品质因数 % beta - 陷波深度参数 (0 <= beta <= 1)。beta=0为最大衰减，beta=1为全通。 % fs - 采样频率 (Hz) % 输出: % b, a - 滤波器差分方程的分子和分母系数向量，a(1)=1。 % % 作者：基于理论推导的MATLAB实现 % 1. 计算预畸变频率 T = 1/fs; K = 2/T; % 双线性变换常数 w0_digital = 2*pi*f0/fs; % 数字角频率 w0_prewarped = K * tan(w0_digital / 2); % 预畸变后的模拟角频率 % 2. 计算中间变量 alpha = w0_prewarped / Q; Omega0_sq = w0_prewarped^2; K_sq = K^2; % 3. 计算公共分母 D D = K_sq + alpha*K + Omega0_sq; % 4. 计算滤波器系数 (已归一化，使得 a(1)=1) b0 = (K_sq + beta*alpha*K + Omega0_sq) / D; b1 = (-2*K_sq + 2*Omega0_sq) / D; b2 = (K_sq - beta*alpha*K + Omega0_sq) / D; a1 = (-2*K_sq + 2*Omega0_sq) / D; % 注意：b1 等于 a1 a2 = (K_sq - alpha*K + Omega0_sq) / D; % 组装系数向量 b = [b0, b1, b2]; a = [1, a1, a2]; % 5. (可选) 绘制频率响应图 if nargout == 0 % 如果无输出参数，则自动绘图 figure('Position', [100, 100, 900, 600]); % 幅频响应 subplot(2,1,1); [h, f] = freqz(b, a, 4096, fs); plot(f, 20*log10(abs(h)), 'LineWidth', 1.5); grid on; xlabel('频率 (Hz)'); ylabel('增益 (dB)'); title(sprintf('三参数陷波滤波器频率响应 (f0=%.1f Hz, Q=%.2f, \\beta=%.2f)', f0, Q, beta)); xlim([0, fs/2]); yline(20*log10(beta), 'r--', 'LineWidth', 1.2); % 标记理论深度线 legend('响应曲线', sprintf('理论深度: %.2f dB', 20*log10(beta)), 'Location', 'best'); % 相频响应 subplot(2,1,2); plot(f, angle(h)*180/pi, 'LineWidth', 1.5); grid on; xlabel('频率 (Hz)'); ylabel('相位 (度)'); xlim([0, fs/2]); end end

代码使用示例与验证：

%% 设计示例 fs = 1000; % 采样率 1kHz f0 = 50; % 陷波频率 50Hz (例如工频干扰) Q = 5; % 中等品质因数 beta = 0.1; % 深度参数，目标在50Hz处衰减到0.1倍 (-20dB) % 获取滤波器系数 [b, a] = three_param_notch(f0, Q, beta, fs); % 生成测试信号：包含10Hz, 50Hz, 100Hz的正弦波 t = 0:1/fs:1-1/fs; % 1秒时长 x = sin(2*pi*10*t) + 0.5*sin(2*pi*50*t) + 0.3*sin(2*pi*100*t); % 使用滤波器处理信号 y = filter(b, a, x); % 绘制结果 figure('Position', [100, 100, 1200, 400]); subplot(1,3,1); plot(t(1:200), x(1:200), 'b'); hold on; plot(t(1:200), y(1:200), 'r', 'LineWidth', 1.5); xlabel('时间 (s)'); ylabel('幅值'); legend('原始信号', '滤波后信号'); title('时域波形对比 (前200个点)'); grid on; subplot(1,3,2); [Pxx, F] = pwelch(x, hanning(256), 128, 256, fs); [Pyy, F] = pwelch(y, hanning(256), 128, 256, fs); plot(F, 10*log10(Pxx), 'b'); hold on; plot(F, 10*log10(Pyy), 'r', 'LineWidth', 1.5); xlabel('频率 (Hz)'); ylabel('功率谱密度 (dB/Hz)'); legend('原始信号谱', '滤波后信号谱'); title('功率谱密度对比'); xlim([0, 150]); grid on; % 标记陷波频率 xline(f0, 'k--', 'LineWidth', 0.8); text(f0+2, -10, sprintf('f0=%dHz', f0), 'FontSize', 9); subplot(1,3,3); % 单独绘制50Hz成分的衰减效果 idx_50hz = find(abs(F-f0) < 2, 1); % 找到50Hz附近的频点 attenuation_db = 10*log10(Pyy(idx_50hz) / Pxx(idx_50hz)); bar(1, attenuation_db, 'FaceColor', [0.8, 0.2, 0.2]); ylabel('衰减 (dB)'); title(sprintf('在%.1fHz处的实际衰减: %.2f dB', F(idx_50hz), attenuation_db)); grid on; ylim([-40, 5]);

运行这段代码，你会看到：

1. 频率响应图上，在50Hz处有一个明显的凹陷，凹陷的底部恰好在我们设定的-20dB（20*log10(0.1)）线附近。
2. 时域波形中，50Hz的成分被显著削弱。
3. 频谱对比图清晰显示50Hz谱峰被抑制。
4. 最后一张图定量显示在50Hz处的衰减接近-20dB，验证了参数 (\beta) 对深度的精确控制。

6. 关键参数影响分析与设计经验

在实际使用三参数陷波滤波器时，理解各个参数对滤波器性能的相互影响至关重要。这能帮助你在调试时快速找到问题所在。

1. 中心频率 (f_0) 与预畸变：这是精度要求最高的参数。务必使用预畸变公式计算。一个快速验证方法是：用freqz函数画出频率响应后，检查-3dB点或谷底点是否在你设定的 (f_0) 上。如果偏差较大，首先检查预畸变计算和采样率 (f_s) 是否正确。对于固定频率的干扰（如50Hz工频），这是一个固定值。对于需要跟踪变化频率的应用（如变速电机谐波），你需要实时更新 (f_0) 并重新计算系数。

2. 品质因数 (Q)：Q值直接决定陷波的“宽度”。其定义是中心频率与-3dB带宽的比值：(Q = f_0 / BW_{-3dB})。

• 高Q值（如Q>10）：陷波非常窄，只滤除极其接近 (f_0) 的频率，对周边频率影响小。适用于抑制一个很纯的单频干扰。但过高的Q值会导致滤波器系数对量化误差非常敏感，在定点DSP或FPGA实现中可能不稳定。
• 低Q值（如Q<2）：陷波很宽，能滤除一个频带内的干扰，但也会损伤该频带内有用的信号。适用于干扰频率有一定波动或带宽的情况。
• 经验选择：通常从 (Q = f_0 / (预期带宽)) 开始估算。例如，要抑制49Hz到51Hz的干扰，带宽约2Hz，则初始Q值可设为50/2=25。然后通过仿真微调。

3. 深度参数 (\beta)：这是三参数滤波器区别于标准陷波器的核心。

• (\beta = 0)：最大衰减。在模拟域理论增益为0，数字域受限于系数精度，通常能达到-60dB到-100dB的衰减。
• (0 < \beta < 1)：可控衰减。例如 (\beta=0.5) 对应-6dB衰减，(\beta=0.1)对应-20dB衰减。这里有一个重要经验：(\beta) 不能太接近1。例如，如果你只想做轻微衰减（如-1dB，对应 (\beta=0.89)），由于双线性变换的非线性以及系数舍入误差，实际频率响应在 (f_0) 处的增益可能与你设定的 (\beta) 有较大出入，且陷波形状可能变得不理想。对于需要浅衰减的场景，建议换用其他结构（如PEAK滤波器或Shelving滤波器）。
• 深度与稳定性的权衡：理论上，只要 (\beta \le 1)，滤波器就是稳定的（因为它是全通和标准陷波的凸组合）。但当 (\beta) 非常小（追求极深陷波）且Q值非常高时，分母系数 (a_2) 会非常接近1，这可能导致在有限精度运算中出现极限环振荡或对输入噪声过于敏感。在FPGA实现时，需要足够的字长来处理这些接近1的系数。

4. 采样频率 (f_s)：

• 奈奎斯特限制：显然，(f_0) 必须小于 (f_s/2)。
• 过采样效应：如果 (f_s) 远高于 (f_0)（例如 (f_s > 20 f_0)），预畸变效应会变得非常微弱，因为 (\tan(\pi f_0/f_s) \approx \pi f_0/f_s)。此时可以近似忽略预畸变，简化计算。但对于音频（fs=44.1kHz）中处理低频（如100Hz以下）干扰，或者电力电子中较低的采样率处理工频，预畸变是必须的。
• 系数更新率：在自适应应用中，如果 (f_0) 变化很快，你需要以多快的频率重新计算并更新系数 (b, a)？这取决于你的系统实时性要求。通常，系数更新率不需要和采样率一样高，可以在一个控制周期内计算好新系数，然后平滑地切换过去，避免输出跳变。

7. 进阶话题：零极点分析与稳定性验证

对于追求深层次理解的工程师，查看滤波器的零极点图能直观判断其特性。三参数陷波滤波器的z域传递函数有两个极点（Poles）和两个零点（Zeros）。

• 极点：由分母多项式 (1 + a_1 z^{-1} + a_2 z^{-2} = 0) 的根决定。它们决定了滤波器的固有频率和衰减速度。对于稳定的滤波器，所有极点必须位于z平面的单位圆内。我们的设计方法（双线性变换）保证了这一点。
• 零点：由分子多项式 (b_0 + b_1 z^{-1} + b_2 z^{-2} = 0) 的根决定。陷波滤波器的零点通常位于单位圆上或单位圆附近，其角度对应着陷波的中心频率 (\theta_0 = 2\pi f_0 / f_s)。

我们可以用MATLAB轻松绘制零极点图并验证稳定性：

%% 零极点分析与稳定性验证 f0 = 50; Q = 5; beta = 0.01; fs = 1000; [b, a] = three_param_notch(f0, Q, beta, fs); figure; zplane(b, a); % 绘制零极点图 title(sprintf('零极点图 (f0=%dHz, Q=%.1f, \\beta=%.3f)', f0, Q, beta)); grid on; % 计算并显示极点模长 poles = roots(a); fprintf('极点位置:\n'); for i = 1:length(poles) fprintf(' 极点 %d: %.6f %+.6fj, 模长 |r| = %.6f\n', i, real(poles(i)), imag(poles(i)), abs(poles(i))); end if all(abs(poles) < 1) fprintf('--> 所有极点均在单位圆内，滤波器稳定。\n'); else fprintf('--> 警告：存在极点位于单位圆上或之外，滤波器不稳定！\n'); end % 计算零点的角度，并反推频率 zeros = roots(b); zero_angles = angle(zeros); % 弧度 zero_freqs = zero_angles * fs / (2*pi); % Hz fprintf('\n零点位置 (对应陷波频率):\n'); for i = 1:length(zeros) if zero_freqs(i) > 0 % 通常取正频率 fprintf(' 零点 %d: 计算频率 ≈ %.2f Hz (理论值 %.1f Hz)\n', i, zero_freqs(i), f0); end end

运行这段代码，你会看到一对共轭零点位于单位圆上（如果 (\beta=0)）或单位圆内（如果 (\beta>0)），其角度对应着陷波频率。极点则位于单位圆内，靠近零点，但更靠近圆心，这保证了滤波器的稳定性。当 (\beta) 趋近于1时，零点向极点靠近，最终重合（全通）。当 (\beta=0) 时，零点在单位圆上，陷波最深。

一个重要的实操心得：在实时系统（尤其是嵌入式系统）中实现时，直接使用上面推导出的差分方程可能会因为系数精度问题（特别是当Q值很高时）导致数值不稳定。一个更鲁棒的做法是使用二阶直接型II（Direct Form II）的转置结构来实现，它对系数量化误差的敏感度通常低于直接型I。在MATLAB中，filter函数内部已经做了优化。但如果自己用C语言在MCU上实现，结构的选择就很重要。