1. 伺服系统中的陷波滤波器基础
在伺服控制系统中，机械谐振和电气噪声是影响系统性能的两大主要因素。作为一名长期从事伺服系统开发的工程师，我经常需要处理这类问题。陷波滤波器（Notch Filter）就像是一个精准的"频率手术刀"，能够在特定频率点产生深度衰减，同时保持其他频段的信号几乎不受影响。
陷波滤波器的核心参数包括：
- 中心频率（f0）：需要抑制的干扰频率点
- 品质因数（Q）：决定滤波器的带宽，Q值越高带宽越窄
- 衰减深度：在中心频率处的衰减程度
在实际工程中，我们通常先在连续时间域设计好滤波器特性，然后通过离散化方法将其转换为数字滤波器。这个过程需要考虑采样率、频率畸变等因素，这也是为什么双线性变换成为工业界的首选方法。
2. 双线性变换原理与实现
2.1 双线性变换的数学基础
双线性变换的本质是将s平面映射到z平面，其核心公式为：
s = (2/T) * (1 - z⁻¹)/(1 + z⁻¹)
这个变换有几个重要特性：
- 将s平面的左半平面映射到z平面的单位圆内，保持稳定性
- 频率响应存在非线性畸变（频率翘曲效应）
- 计算相对简单，适合实时系统实现
我在多个伺服项目中发现，当采样频率至少是目标频率的10倍时，双线性变换的效果最佳。例如对于50Hz的陷波频率，采样率不应低于500Hz。
2.2 Python实现示例
下面是一个更完整的Python实现，包含了频率响应绘制：
python复制import numpy as np
import control
import matplotlib.pyplot as plt
# 连续时间陷波滤波器参数
f0 = 50 # Hz
Q = 15
omega0 = 2 * np.pi * f0
# 采样参数
fs = 1000 # 采样频率1kHz
T = 1/fs
# 连续时间传递函数
num = [1, 0, omega0**2]
den = [1, omega0/Q, omega0**2]
sys_c = control.TransferFunction(num, den)
# 双线性变换离散化
sys_d = control.sample_system(sys_c, T, method='bilinear')
# 频率响应分析
w = np.logspace(1, 3, 1000) # 10Hz到1kHz
mag_c, phase_c, omega_c = control.bode(sys_c, w, plot=False)
mag_d, phase_d, omega_d = control.bode(sys_d, w, plot=False)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.semilogx(omega_c/(2*np.pi), 20*np.log10(mag_c), label='Continuous')
plt.semilogx(omega_d/(2*np.pi), 20*np.log10(mag_d), label='Discrete')
plt.axvline(f0, color='r', linestyle='--', alpha=0.5)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude [dB]')
plt.grid(which='both')
plt.legend()
plt.title('Notch Filter Frequency Response')
plt.show()
注意：实际工程中，建议使用scipy.signal中的iirnotch函数，它已经内置了频率补偿机制。
3. 频率补偿技术详解
3.1 频率畸变现象分析
双线性变换导致的频率畸变可以用以下公式描述：
ω_d = (2/T) * tan(ω_c * T/2)
这意味着：
- 低频段畸变较小
- 接近Nyquist频率时畸变严重
- 实际陷波频率会低于设计频率
我在一个电机控制项目中实测发现，设计50Hz的陷波器，实际陷波点在47Hz左右，这会导致抑制效果大打折扣。
3.2 实用补偿方法
经过多次实验，我总结出一个简单有效的补偿公式：
f_adjusted = f_design / cos(π * f_design / fs)
对应的Python实现：
python复制def freq_prewarp(f_design, fs):
return f_design / np.cos(np.pi * f_design / fs)
f0_comp = freq_prewarp(f0, fs)
omega0_comp = 2 * np.pi * f0_comp
# 使用补偿频率重新设计
num_comp = [1, 0, omega0_comp**2]
den_comp = [1, omega0_comp/Q, omega0_comp**2]
sys_c_comp = control.TransferFunction(num_comp, den_comp)
sys_d_comp = control.sample_system(sys_c_comp, T, method='bilinear')
这种方法不需要复杂的计算，在大多数工业应用中已经足够精确。我的经验是，当f_design < fs/5时，补偿误差可以控制在1%以内。
4. 工程实现与调试技巧
4.1 离散化实现形式
离散化后的传递函数通常表示为：
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]
在实际DSP实现时，我推荐使用直接II型结构，它只需要两个延迟单元，计算效率最高。
4.2 参数调试经验
-
Q值选择：
- 机械谐振：Q=10~30
- 电气噪声：Q=5~15
- 太高会导致相位突变，影响稳定性
-
采样率选择：
- 最低要求：fs > 4*f0
- 推荐：fs > 10*f0
- 太高会增加计算负担
-
衰减深度调整：
可以通过在分子项添加增益因子实现：
H(z) = K * (b0 + b1 z⁻¹ + b2 z⁻²) / (1 + a1 z⁻¹ + a2 z⁻²)
K=0.5可实现-6dB衰减
5. 常见问题与解决方案
5.1 陷波频率偏移
现象：实际抑制频率与设计值不符
可能原因：
- 频率补偿不准确
- 采样率过低
解决方案： - 重新校准补偿系数
- 提高采样率
- 使用扫频信号实测响应
5.2 滤波器不稳定
现象：输出发散或振荡
可能原因：
- 量化误差累积
- 系数超出范围
解决方案： - 使用更高精度算法（32位浮点）
- 检查系数是否满足|a1|+|a2|<1的稳定性条件
5.3 相位失真严重
现象：系统相位裕度降低
可能原因：
- Q值过高
- 陷波频率接近截止频率
解决方案： - 降低Q值
- 调整陷波频率位置
- 考虑使用多个低Q值陷波器串联
6. 进阶应用：自适应陷波滤波器
在转速变化的场合，固定频率的陷波器效果有限。我最近在一个风电项目中实现了自适应陷波器，核心思路是：
- 实时检测干扰频率（如通过FFT或PLL）
- 动态更新滤波器系数
- 平滑过渡避免突变
关键代码片段：
python复制def update_notch_params(f_new, Q, fs):
f_comp = freq_prewarp(f_new, fs)
omega = 2 * np.pi * f_comp
b = [1, 0, omega**2]
a = [1, omega/Q, omega**2]
return bilinear_transform(b, a, fs)
# 在主循环中调用
while True:
f_est = estimate_frequency(input_signal)
b, a = update_notch_params(f_est, Q=10, fs=1000)
output_signal = apply_filter(b, a, input_signal)
这种实现方式在变速工况下表现优异，但需要注意更新速率不宜过快，通常控制在几个毫秒量级。