陷波滤波器（Notch Filter）的离散化设计
符号说明
w b w w_{bw} wbw - 陷波宽度，单位： r a d / s rad/s rad/s
w c w_c wc - 陷波中央频率，单位： r a d / s rad/s rad/s
T s T_s Ts - 离散化采样时间，单位： s s s
y ( k ) y(k) y(k) - 第k次输出信号
r ( k ) r(k) r(k) - 第k次输入信号
概述
本文以如下二阶陷波滤波器的传递函数为例，分别简要介绍双线性变换和零极点匹配方法的离散化：
G
(
s
)
=
s
2
+
w
c
2
s
2
+
w
b
w
s
+
w
c
2
(0-1)
G(s) = \frac{s^2 + w_c^2}{s^2 + w_{bw}s + w_c^2} \tag{0-1}
G(s)=s2+wbws+wc2s2+wc2(0-1)
双线性变换（Tustin）方法
双线性变换本质是一种数值积分法，采用梯形方法来近似计算积分。经过简要推导可以得到：
s
=
2
T
s
z
−
1
z
+
1
(1-1)
s = \frac{2}{T_s} \frac{z-1}{z+1} \tag{1-1}
s=Ts2z+1z−1(1-1)
即使用式(1-1)代入式(0-1)即可得到陷波滤波器的离散化z域方程：
G
(
z
)
=
4
T
s
(
z
−
1
)
2
(
z
+
1
)
2
+
w
c
2
4
T
s
(
z
−
1
)
2
(
z
+
1
)
2
+
w
b
w
2
T
s
z
−
1
z
+
1
+
w
c
2
=
4
+
w
c
2
T
s
2
⏞
a
0
+
(
2
T
s
2
w
c
2
−
8
)
⏞
a
1
z
−
1
+
(
4
+
w
c
2
T
s
2
)
⏞
a
2
z
−
2
4
+
w
c
2
T
s
2
+
2
w
b
w
T
s
⏟
b
0
+
(
2
T
s
2
w
c
2
−
8
)
⏟
b
1
z
−
1
+
(
4
+
w
c
2
T
s
2
−
2
w
b
w
T
s
)
⏟
b
2
z
−
2
(1-2)
G(z) = \frac{\frac{4}{Ts}\frac{(z-1)^2}{(z+1)^2}+w_c^2}{\frac{4}{Ts}\frac{(z-1)^2}{(z+1)^2}+w_{bw}\frac{2}{Ts}\frac{z-1}{z+1}+w_c^2} \\ = \frac{\overbrace{4+w_c^2T_s^2}^{a_0} + \overbrace{(2T_s^2w_c^2-8)}^{a_1}z^{-1} + \overbrace{(4+w_c^2T_s^2)}^{a_2}z^{-2}}{\underbrace{4+w_c^2T_s^2+2w_{bw}T_s}_{b_0} + \underbrace{(2T_s^2w_c^2-8)}_{b_1}z^{-1} + \underbrace{(4+w_c^2T_s^2-2w_{bw}T_s)}_{b_2}z^{-2}} \tag{1-2}
G(z)=Ts4(z+1)2(z−1)2+wbwTs2z+1z−1+wc2Ts4(z+1)2(z−1)2+wc2=b0
4+wc2Ts2+2wbwTs+b1
(2Ts2wc2−8)z−1+b2
(4+wc2Ts2−2wbwTs)z−24+wc2Ts2
a0+(2Ts2wc2−8)
a1z−1+(4+wc2Ts2)
a2z−2(1-2)
对式(1-2)作变量替代可得到：
G
(
z
)
=
a
0
+
a
1
z
−
1
+
a
2
z
−
2
b
0
+
b
1
z
−
1
+
b
2
z
−
2
(1-3)
G(z) = \frac{a_0 + a_1z^{-1} + a_2z^{-2}}{b_0 + b_1z^{-1} + b_2z^{-2}} \tag{1-3}
G(z)=b0+b1z−1+b2z−2a0+a1z−1+a2z−2(1-3)
其中：
- a 0 = 4 + w c 2 T s 2 a_0 = 4+w_c^2T_s^2 a0=4+wc2Ts2
- a 1 = 2 T s 2 w c 2 − 8 a_1 = 2T_s^2w_c^2-8 a1=2Ts2wc2−8
- a 2 = 4 + w c 2 T s 2 = a 0 a_2 = 4+w_c^2T_s^2 = a_0 a2=4+wc2Ts2=a0
- b 0 = 4 + w c 2 T s 2 + 2 w b w T s b_0 = 4+w_c^2T_s^2+2w_{bw}T_s b0=4+wc2Ts2+2wbwTs
- b 1 = 2 T s 2 w c 2 − 8 = a 1 b_1 = 2T_s^2w_c^2-8 = a_1 b1=2Ts2wc2−8=a1
- b 2 = 4 + w c 2 T s 2 − 2 w b w T s b_2 = 4+w_c^2T_s^2-2w_{bw}T_s b2=4+wc2Ts2−2wbwTs
由式(1-3)可得出陷波滤波器基于Tustin方法的离散化差分方程为：
y
(
k
)
=
a
0
b
0
r
(
k
)
+
a
1
b
0
r
(
k
−
1
)
+
a
2
b
0
r
(
k
−
2
)
−
b
1
b
0
y
(
k
−
1
)
−
b
2
b
0
y
(
k
−
2
)
(1-4)
y(k) = \frac{a_0}{b_0}r(k)+\frac{a_1}{b_0}r(k-1)+\frac{a_2}{b_0}r(k-2)-\frac{b1}{b0}y(k-1)-\frac{b2}{b0}y(k-2) \tag{1-4}
y(k)=b0a0r(k)+b0a1r(k−1)+b0a2r(k−2)−b0b1y(k−1)−b0b2y(k−2)(1-4)
值得注意的是，若使用后向欧拉法离散，只需要将式(1-1)更换成下式即可：
s
=
1
−
z
−
1
T
s
(1-5)
s = \frac{1-z^{-1}}{T_s} \tag{1-5}
s=Ts1−z−1(1-5)
零极点匹配方法
零极点匹配是指将s域中的零极点一一对应到z域的零极点上，并计算出增益即可。零极点对应关系由下式给出：
z
i
=
e
s
i
T
s
(2-1)
z_i = e^{s_iT_s} \tag{2-1}
zi=esiTs(2-1)
其中下标i表示第i个零点或极点。可将式(0-1)表示成：
G
(
s
)
=
K
s
(
s
−
j
w
c
)
(
s
+
j
w
c
)
(
s
−
p
1
)
(
s
−
p
2
)
(2-2)
G(s) = K_s\frac{(s-jw_c)(s+jw_c)}{(s-p_1)(s-p_2)} \tag{2-2}
G(s)=Ks(s−p1)(s−p2)(s−jwc)(s+jwc)(2-2)
其中：
- K s = 1 K_s = 1 Ks=1
- p 1 , 2 = α ± β p_{1,2} = \alpha \pm \beta p1,2=α±β，为陷波滤波器的极点
- α = − w b w 2 , β = 4 w c 2 − w b w 2 2 \alpha = -\frac{w_{bw}}{2}, \beta = \frac{\sqrt{4w_c^2-w_{bw}^2}}{2} α=−2wbw,β=24wc2−wbw2
将式(2-2)中的零极点利用式(2-1)替换，借助欧拉公式，可得到z域下的方程如下：
G
(
z
)
=
K
z
(
z
−
e
j
w
c
T
s
)
(
z
−
e
−
j
w
c
T
s
)
(
z
−
e
p
1
T
s
)
(
z
−
e
p
2
T
s
)
=
K
z
(
z
−
e
j
w
c
T
s
)
(
z
−
e
−
j
w
c
T
s
)
(
z
−
r
e
j
β
T
s
)
(
z
−
r
e
−
j
β
T
s
)
=
K
z
z
2
−
2
c
o
s
(
w
c
T
s
)
z
+
1
z
2
−
2
r
c
o
s
(
β
T
s
)
z
+
r
2
(2-3)
G(z) = K_z\frac{(z-e^{jw_cT_s})(z-e^{-jw_cT_s})}{(z-e^{p_1T_s})(z-e^{p_2T_s})} \\ = K_z\frac{(z-e^{jw_cT_s})(z-e^{-jw_cT_s})}{(z-re^{j\beta T_s})(z-re^{-j\beta T_s})} \\ = K_z\frac{z^2 - 2cos(w_cT_s)z + 1}{z^2 - 2rcos(\beta T_s)z + r^2} \tag{2-3}
G(z)=Kz(z−ep1Ts)(z−ep2Ts)(z−ejwcTs)(z−e−jwcTs)=Kz(z−rejβTs)(z−re−jβTs)(z−ejwcTs)(z−e−jwcTs)=Kzz2−2rcos(βTs)z+r2z2−2cos(wcTs)z+1(2-3)
其中：
- r = e α T s r = e^{\alpha T_s} r=eαTs
当s域中为0时，z域对应为1，因此通过式(2-2)与式(2-3)可得到下式：
G
(
s
)
∣
s
=
0
=
G
(
z
)
∣
z
=
1
⇒
1
=
K
z
2
−
2
c
o
s
(
w
c
T
s
)
1
−
2
r
c
o
s
(
β
T
s
)
+
r
2
⇒
K
z
=
1
−
2
r
c
o
s
(
β
T
s
)
+
r
2
2
−
2
c
o
s
(
w
c
T
s
)
(2-4)
G(s)\mid_{s=0} = G(z)\mid_{z=1} \\ \Rightarrow 1 = K_z\frac{2 - 2cos(w_cT_s)}{1-2rcos(\beta T_s) + r^2} \\ \Rightarrow K_z = \frac{1-2rcos(\beta T_s) + r^2}{2-2cos(w_cT_s)} \tag{2-4}
G(s)∣s=0=G(z)∣z=1⇒1=Kz1−2rcos(βTs)+r22−2cos(wcTs)⇒Kz=2−2cos(wcTs)1−2rcos(βTs)+r2(2-4)
使用变量替代，可将式(2-3)化简得到：
G
(
z
)
=
K
z
⏞
a
0
z
2
−
2
K
z
c
o
s
(
w
c
T
s
)
⏞
a
1
z
+
K
z
⏞
a
2
z
2
−
2
r
c
o
s
(
β
T
s
)
⏟
b
1
z
+
r
2
⏟
−
b
2
=
a
0
z
2
+
a
1
z
+
a
2
z
2
−
b
1
−
b
2
(2-5)
G(z) = \frac{\overbrace{K_z}^{a_0}z^2\overbrace{ - 2K_zcos(w_cT_s)}^{a_1}z + \overbrace{K_z}^{a_2}}{z^2 - \underbrace{2rcos(\beta T_s)}_{b_1}z + \underbrace{r^2}_{-b_2}} \\ = \frac{a_0z^2 + a_1z + a_2}{z^2 - b_1 - b_2} \tag{2-5}
G(z)=z2−b1
2rcos(βTs)z+−b2
r2Kz
a0z2−2Kzcos(wcTs)
a1z+Kz
a2=z2−b1−b2a0z2+a1z+a2(2-5)
其中：
- a 0 = K z a_0 = K_z a0=Kz
- a 1 = − 2 K z c o s ( w c T s ) a_1 = -2K_zcos(w_cT_s) a1=−2Kzcos(wcTs)
- a 2 = K z = a 0 a_2 = K_z = a_0 a2=Kz=a0
- b 1 = 2 r c o s ( β T s ) b_1 = 2rcos(\beta T_s) b1=2rcos(βTs)
- b 2 = − r 2 b_2 = -r^2 b2=−r2
由式(2-5)可得出可得出陷波滤波器基于零极点匹配方法的离散化差分方程为：
y
(
k
)
=
a
0
r
(
k
)
+
a
1
r
(
k
−
1
)
+
a
2
r
(
k
−
2
)
+
b
1
y
(
k
−
1
)
+
b
2
y
(
k
−
2
)
(2-6)
y(k) = a_0r(k)+a_1r(k-1)+a_2r(k-2)+b_1y(k-1)+b_2y(k-2) \tag{2-6}
y(k)=a0r(k)+a1r(k−1)+a2r(k−2)+b1y(k−1)+b2y(k−2)(2-6)
参考资料
- 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');
628
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
