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)。(\epsilon0) 代表完全抑制(\epsilon1) 则相当于全通无任何滤波效果。一种常见且数学上优雅的引入方式是将分子项进行加权[ 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)让我们来分析一下当 (\beta 1) 时分子和分母完全相同(H_{3p}(s) 1)滤波器变成一个全通网络没有任何滤波效果。这对应深度最浅无衰减的情况。当 (\beta 0) 时分子变回 (s^2 \omega_0^2)这就是标准的二阶陷波滤波器在中心频率处衰减最大。这对应深度最深的情况。当 (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 Warpings域的模拟频率 (\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) ]这是整个离散化过程中最关键、也最容易出错的一步。很多初学者实现的滤波器频率不准问题就出在忽略了预畸变。我们的离散化步骤如下确定数字指标给定采样频率 (f_s)期望陷波的数字中心频率 (f_d f_0)品质因数 (Q)陷波深度 (\beta)。预畸变计算计算预畸变后的模拟角频率 (\omega_0 2 \pi f_0 \frac{2}{T} \tan(\pi f_0 / f_s))。构造模拟传递函数将 (\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} ]应用双线性变换将 (s \frac{2}{T} \cdot \frac{1 - z^{-1}}{1 z^{-1}}) 代入上式。整理为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_01) [ \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_01, b_1a_1, b_2a_2)滤波器确为全通当 (\beta0) 时(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)。beta0为最大衰减beta1为全通。 % 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(f02, -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]);运行这段代码你会看到频率响应图上在50Hz处有一个明显的凹陷凹陷的底部恰好在我们设定的-20dB20*log10(0.1)线附近。时域波形中50Hz的成分被显著削弱。频谱对比图清晰显示50Hz谱峰被抑制。最后一张图定量显示在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值如Q10陷波非常窄只滤除极其接近 (f_0) 的频率对周边频率影响小。适用于抑制一个很纯的单频干扰。但过高的Q值会导致滤波器系数对量化误差非常敏感在定点DSP或FPGA实现中可能不稳定。低Q值如Q2陷波很宽能滤除一个频带内的干扰但也会损伤该频带内有用的信号。适用于干扰频率有一定波动或带宽的情况。经验选择通常从 (Q f_0 / (预期带宽)) 开始估算。例如要抑制49Hz到51Hz的干扰带宽约2Hz则初始Q值可设为50/225。然后通过仿真微调。3. 深度参数 (\beta)这是三参数滤波器区别于标准陷波器的核心。(\beta 0)最大衰减。在模拟域理论增益为0数字域受限于系数精度通常能达到-60dB到-100dB的衰减。(0 \beta 1)可控衰减。例如 (\beta0.5) 对应-6dB衰减(\beta0.1)对应-20dB衰减。这里有一个重要经验(\beta) 不能太接近1。例如如果你只想做轻微衰减如-1dB对应 (\beta0.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)。此时可以近似忽略预畸变简化计算。但对于音频fs44.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运行这段代码你会看到一对共轭零点位于单位圆上如果 (\beta0)或单位圆内如果 (\beta0)其角度对应着陷波频率。极点则位于单位圆内靠近零点但更靠近圆心这保证了滤波器的稳定性。当 (\beta) 趋近于1时零点向极点靠近最终重合全通。当 (\beta0) 时零点在单位圆上陷波最深。一个重要的实操心得在实时系统尤其是嵌入式系统中实现时直接使用上面推导出的差分方程可能会因为系数精度问题特别是当Q值很高时导致数值不稳定。一个更鲁棒的做法是使用二阶直接型IIDirect Form II的转置结构来实现它对系数量化误差的敏感度通常低于直接型I。在MATLAB中filter函数内部已经做了优化。但如果自己用C语言在MCU上实现结构的选择就很重要。
