陷波滤波器（Notch Filter）的离散化设计

作者: SmartFish | 来源:发表于2021-05-13 11:52 被阅读0次

陷波滤波器（Notch Filter）的离散化设计

符号说明

$w_{bw}$ - 陷波宽度，单位： $rad/s$

$w_c$ - 陷波中央频率，单位： $rad/s$

$T_s$ - 离散化采样时间，单位： $s$

$y(k)$ - 第k次输出信号

$r(k)$ - 第k次输入信号

概述

本文以如下二阶陷波滤波器的传递函数为例，分别简要介绍双线性变换和零极点匹配方法的离散化：
$G(s) = \frac{s^2 + w_c^2}{s^2 + w_{bw}s + w_c^2} \tag{0-1}$

双线性变换（Tustin）方法

双线性变换本质是一种数值积分法，采用梯形方法来近似计算积分。经过简要推导可以得到：
$s = \frac{2}{T_s} \frac{z-1}{z+1} \tag{1-1}$
即使用式(1-1)代入式(0-1)即可得到陷波滤波器的离散化z域方程：
$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}$
对式(1-2)作变量替代可得到：
$G(z) = \frac{a_0 + a_1z^{-1} + a_2z^{-2}}{b_0 + b_1z^{-1} + b_2z^{-2}} \tag{1-3}$
其中：

$a_0 = 4+w_c^2T_s^2$
$a_1 = 2T_s^2w_c^2-8$
$a_2 = 4+w_c^2T_s^2 = a_0$
$b_0 = 4+w_c^2T_s^2+2w_{bw}T_s$
$b_1 = 2T_s^2w_c^2-8 = a_1$
$b_2 = 4+w_c^2T_s^2-2w_{bw}T_s$

由式(1-3)可得出陷波滤波器基于Tustin方法的离散化差分方程为：
$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}$
值得注意的是，若使用后向欧拉法离散，只需要将式(1-1)更换成下式即可：
$s = \frac{1-z^{-1}}{T_s} \tag{1-5}$

零极点匹配方法

零极点匹配是指将s域中的零极点一一对应到z域的零极点上，并计算出增益即可。零极点对应关系由下式给出：
$z_i = e^{s_iT_s} \tag{2-1}$
其中下标i表示第i个零点或极点。可将式(0-1)表示成：
$G(s) = K_s\frac{(s-jw_c)(s+jw_c)}{(s-p_1)(s-p_2)} \tag{2-2}$
其中：

$K_s = 1$
$p_{1,2} = \alpha \pm \beta$ ，为陷波滤波器的极点
$\alpha = -\frac{w_{bw}}{2}, \beta = \frac{\sqrt{4w_c^2-w_{bw}^2}}{2}$

将式(2-2)中的零极点利用式(2-1)替换，借助欧拉公式，可得到z域下的方程如下：
$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}$
其中：

$r = e^{\alpha T_s}$

当s域中为0时，z域对应为1，因此通过式(2-2)与式(2-3)可得到下式：
$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}$
使用变量替代，可将式(2-3)化简得到：
$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}$
其中：

$a_0 = K_z$
$a_1 = -2K_zcos(w_cT_s)$
$a_2 = K_z = a_0$
$b_1 = 2rcos(\beta T_s)$
$b_2 = -r^2$

由式(2-5)可得出可得出陷波滤波器基于零极点匹配方法的离散化差分方程为：
$y(k) = a_0r(k)+a_1r(k-1)+a_2r(k-2)+b_1y(k-1)+b_2y(k-2) \tag{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');

陷波滤波器（Notch Filter）的离散化设计