陷波滤波器（Notch Filter）的离散化设计
符号说明
- 陷波宽度，单位：
- 陷波中央频率，单位：
- 离散化采样时间，单位：
- 第k次输出信号
- 第k次输入信号
概述
本文以如下二阶陷波滤波器的传递函数为例，分别简要介绍双线性变换和零极点匹配方法的离散化：
双线性变换（Tustin）方法
双线性变换本质是一种数值积分法，采用梯形方法来近似计算积分。经过简要推导可以得到：
即使用式(1-1)代入式(0-1)即可得到陷波滤波器的离散化z域方程：
对式(1-2)作变量替代可得到：
其中：
由式(1-3)可得出陷波滤波器基于Tustin方法的离散化差分方程为：
值得注意的是，若使用后向欧拉法离散，只需要将式(1-1)更换成下式即可：
零极点匹配方法
零极点匹配是指将s域中的零极点一一对应到z域的零极点上，并计算出增益即可。零极点对应关系由下式给出：
其中下标i表示第i个零点或极点。可将式(0-1)表示成：
其中：
-
，为陷波滤波器的极点
将式(2-2)中的零极点利用式(2-1)替换，借助欧拉公式，可得到z域下的方程如下：
其中：
当s域中为0时，z域对应为1，因此通过式(2-2)与式(2-3)可得到下式：
使用变量替代，可将式(2-3)化简得到：
其中：
由式(2-5)可得出可得出陷波滤波器基于零极点匹配方法的离散化差分方程为：
参考资料
- Control Systems in Practice, Part 5: A Better Way to Think About a Notch Filter
- 刘建昌、关守平、周玮主编.《计算机控制系统》（第二版）.北京：科学出版社. 2016.8
附录
附两种方法的matlab测试源码
clear all;
close all;
clc;
fc = 100;
fbw = 40;
wc = 2 * pi * fc;
wbw = 2 * pi * fbw;
Ts = 0.001;
a = [1 0 wc^2];
b = [1 wbw wc^2];
sys = tf(a, b);
sysd_tustin = c2d(sys, Ts, 'tustin');
sysd_matched = c2d(sys, Ts, 'matched');
%% tustin test
a0 = 4 + 4 * pi^2 * fc^2 * Ts^2;
a1 = 8 * pi^2 * fc^2 * Ts^2 - 8;
a2 = a0;
b0 = a0 + 4 * pi * fbw * Ts;
b1 = a1;
b2 = a0 - 4 * pi * fbw * Ts;
ad = [a0 a1 a2]./b0;
bd = [b0 b1 b2]./b0;
sysd_tustin_test = tf(ad, bd, Ts);
%% zpm test
alpha = -wbw / 2;
beta = sqrt(4*wc^2-wbw^2) / 2;
r = exp(alpha * Ts);
Kz = (1 - 2 * r * cos(beta * Ts) + r^2) / (2 - 2 * cos(wc * Ts));
a0 = Kz;
a1 = -2 * Kz * cos(wc * Ts);
a2 = Kz;
b0 = 1;
b1 = 2 * r * cos(beta * Ts);
b2 = -r^2;
ad = [a0 a1 a2];
bd = [b0 -b1 -b2];
sysd_matched_test = tf(ad, bd, Ts);
%% figure
figure(1);
P=bodeoptions;
P.FreqUnits = 'Hz';
bode(sys, P);
grid on;
title('sys');
figure(2);
bode(sysd_tustin_test, P);
grid on;
title('sysd\_tustin\_test');
figure(3);
bode(sysd_matched_test, P);
grid on;
title('sysd\_matched\_test');
%% signal test
t = 0:Ts:1;
f0 = 10;
f1 = 0;
f2 = fc;
f3 = 0;
r = sin(2*pi*f0*t) + sin(2*pi*f1*t) + sin(2*pi*f2*t) + sin(2*pi*f3*t);
y_sys_filtered = lsim(sys, r, t);
y_sysd_tustin_filtered = dlsim(sysd_tustin_test.num, sysd_tustin_test.den, r);
y_sysd_matched_filtered = dlsim(sysd_matched_test.num, sysd_matched_test.den, r);
%% display
figure(4);
lw = 2;
plot(t, r);
hold on;
grid on;
plot(t, y_sys_filtered, 'LineWidth',lw);
hold on;
plot(t, y_sysd_tustin_filtered, 'LineWidth',lw);
hold on;
plot(t, y_sysd_matched_filtered, 'LineWidth',lw);
hold on;
legend('Input', 'y\_sys\_filtered', 'y\_sysd\_tustin\_filtered', 'y\_sysd\_matched\_filtered');
title('filtered signal');
