<div id="content_views" class="htmledit_views">
<h2><a name="t0"></a></h2>
一、陷波器在连续域的传递函数
1、最基本的陷波器传函
(1)
其中，wo是所谓“中心频率”，也就是你想要“陷掉”的频率。而 ζ 则是“陷阱”的宽度。
根据公式可以发现，当输入信号频率很小（s=0）或者很大（ s=+∞）的时候，上面式子的值是1；当输入信号频率刚好等于 s=jωo的时候，分子是0，所以增益变成0，那这个频率的信号当然就全都被衰减掉了。
由上图可见，ζ越大，则弦波带宽越宽，但弦波频率处的衰减越小。
2、三参数陷波器传函
(2)
其中，ωo是陷波频率（即凹陷的中心频率），ζ1和ζ2是陷波系数
陷波滤波器重点关注的参数一般有三个：
（1）陷波频率（ωo rad/s可转换Hz）
（2）陷波深度（depth为衰减倍数）
例如对于100Hz频率处的衰减深度是100，那经过该滤波器后，幅值衰减100倍。
（3）陷波宽度（△f单位Hz）
即中心频率两侧，幅值衰减-3dB时，对应的两个频率的差值。
ζ1和ζ2与陷波深度depth和陷波宽度△f（Hz）的关系表示如下：
二、陷波器传函在MATLAB中的表达及其离散化
1、以最基本的陷波器传函为例
a、MATLAB中编写如下m文件：
-
syms w0 s Ts z zeta % 定义符号变量
-
G1 =(s^
2+w0^
2)/(s^
2+
2*w0*zeta*s+w0^
2) %传递函数
-
sys_s2c =
2*(z
-1)/Ts/(z+
1);
-
G2 =
subs(G1,s,sys_s2c) %离散化 tustin变换
-
G3 =
collect(G2,z) % 将表达式G2中的以z为变量合并相同次幂；
得到连续域和Z域的传递函数如下：
G1 =(s^2 + w0^2)/(s^2 + 2*zeta*s*w0 + w0^2)
G2 =(w0^2 + (2*z - 2)^2/(Ts^2*(z + 1)^2))/(w0^2 + (2*z - 2)^2/(Ts^2*(z + 1)^2) + (2*w0*zeta*(2*z - 2))/(Ts*(z + 1)))
G3 =((Ts^2*w0^2 + 4)*z^2 + (2*Ts^2*w0^2 - 8)*z + Ts^2*w0^2 + 4)/((Ts^2*w0^2 + 4*zeta*Ts*w0 + 4)*z^2 + (2*Ts^2*w0^2 - 8)*z + Ts^2*w0^2 - 4*zeta*Ts*w0 + 4)
根据离散化的方法：
分子分母同时除以A0，得到差分方程系数
差分方程为：
y(k) = b0*x(k) + b1*x(k-1) +b2*x(k-2) - a1*y(k-1) - a2*y(k-2);
则
B0 = Ts^2*w0^2 + 4;
B1 = 2*Ts^2*w0^2 - 8;
B2 = B0; （3）
A0 = Ts^2*w0^2 + 4*zeta*Ts*w0 + 4;
A1 = B1;
A2 = Ts^2*w0^2 - 4*zeta*Ts*w0 + 4;
可得差分方程系数
a0 = 1
b0 = B0/A0
b1 = B1/A0 （4）
b2 = B2/A0
a1 = A1/A0
a2 = A2/A0
b、代入实际参数，求得差分方程系数
-
f0 =
100; %目标频率
-
w0 =
2*pi*f0;
-
Ts =
10e-6; %离散化周期
-
zeta =
0.5; %带宽
-
-
% Control object funciton
-
gxnum2 = [
1
0 w0^
2]; % 传函分子
-
gxden3 = [
1
2*w0*zeta w0^
2]; % 传函分母
-
sys2 =
tf(gxnum2, gxden3) % 传函
-
zpk(sys2,
's'); % 将传函用零极点格式表示
-
dsys2 =
c2d(sys2, Ts,
'tustin') % s到z，即，连续到离散域的变换
-
zpk(dsys2,
'z'); % 将传函用零极点格式表示
-
[num2, den2] =
tfdata(dsys2,
'v'); % 提取传递函数分子、分母的z次幂的系数，并保存到数组
得到传递函数如下：
需要注意的是，在上图Z域传递函数dsys2的各项系数是经四舍五入的，如果直接用这些系数到Simulink或编写程序进行仿真，得到的结果是错误的！！！需要使用更加精确的系数。这些系数可以将式3和式4代入MATLAB中直接求解得到。
b0 = 0.996868276853708 b1 = -1.993697199313698
b2 = 0.996868276853708 a1 = -1.993697199313698
a2 = 0.993736553707416
2、以三参数陷波器为例
三、仿真验证
1、使用经四舍五入的系数仿真
结果错误
2、使用正确的系数
结果正确：
四、C程序验证
1、方式一
直接将差分方程用代码呈现：
-
double ADValue;
-
-
//Notching Filter Coefficient
-
#define Notching_filter_100Hz_a0 1
-
#define Notching_filter_100Hz_a1 -1.998704711158930
-
#define Notching_filter_100Hz_a2 0.998744164397945
-
#define Notching_filter_100Hz_b0 0.999372082198973
-
#define Notching_filter_100Hz_b1 -1.998704711158930
-
#define Notching_filter_100Hz_b2 0.999372082198973
-
-
// Control law 2p2z data define
-
typedef
struct
IIR_2OR_DATA_TAG{
-
double coeff_a0;
-
double coeff_a1;
-
double coeff_a2;
-
double coeff_b0;
-
double coeff_b1;
-
double coeff_b2;
-
-
double filter_out;
-
double filter_y1;
-
double filter_y2;
-
double filter_u1;
-
double filter_u2;
-
} IIR_2OR_DATA_DEF;
-
-
IIR_2OR_DATA_DEF notching_data = {
-
Notching_filter_100Hz_a0,
-
Notching_filter_100Hz_a1,
-
Notching_filter_100Hz_a2,
-
Notching_filter_100Hz_b0,
-
Notching_filter_100Hz_b1,
-
Notching_filter_100Hz_b2,
-
0,
0,
0,
0 ,
0
-
};
-
-
double iir_2or_func(IIR_2OR_DATA_DEF *filter_date, double target)
-
{
-
//
-
//y(out) = b0*x(k) + b1*x(k-1) +b2*x(k-2) - a1*y(k-1) - a2*y(k-2);
-
//
-
filter_date->filter_out = (filter_date->coeff_b0 * (target)) \
-
+ (filter_date->coeff_b1 * filter_date->filter_u1) \
-
+ (filter_date->coeff_b2 * filter_date->filter_u2) \
-
- (filter_date->coeff_a1 * filter_date->filter_y1) \
-
- (filter_date->coeff_a2 * filter_date->filter_y2);
-
-
//
-
// Update last data
-
//
-
filter_date->filter_y2 = filter_date->filter_y1;
-
filter_date->filter_y1 = filter_date->filter_out;
-
filter_date->filter_u2 = filter_date->filter_u1;
-
filter_date->filter_u1 = target;
-
-
//
-
// Return Value
-
//
-
return(filter_date->filter_out);
-
}
-
得到结果正确：
2、方式二
参考文献3中的方法，增加一个中间变量w，来实现。
-
double ADValue;
-
-
//Notching Filter Coefficient
-
-
#define Notching_filter_100Hz_GAIN 1
-
#define Notching_filter_100Hz_A1 -1.998704711158930
-
#define Notching_filter_100Hz_A2 0.998744164397945
-
#define Notching_filter_100Hz_B0 0.999372082198973
-
#define Notching_filter_100Hz_B1 -1.998704711158930
-
#define Notching_filter_100Hz_B2 0.999372082198973
-
-
-
// Control law 2p2z data define
-
typedef
struct
IIR_2OR_DATA_TAG{
-
double coeff_GAIN;
-
double coeff_B0;
-
double coeff_B1;
-
double coeff_B2;
-
double coeff_A1;
-
double coeff_A2;
-
-
double filter_out;
-
double filter_W0;
-
double filter_W1;
-
double filter_W2;
-
} IIR_2OR_DATA_DEF;
-
-
IIR_2OR_DATA_DEF notching_data = {
-
Notching_filter_100Hz_GAIN,
-
Notching_filter_100Hz_B0,
-
Notching_filter_100Hz_B1,
-
Notching_filter_100Hz_B2,
-
Notching_filter_100Hz_A1,
-
Notching_filter_100Hz_A2,
-
0,
0,
0,
0
-
};
-
-
double iir_2or_func(IIR_2OR_DATA_DEF *filter_date, double target)
-
{
-
//
-
// w0 = x(0) - A1 * W1 - A2 * W2
-
//
-
filter_date->filter_W0 = (target) - filter_date->coeff_A1 * filter_date->filter_W1 \
-
- filter_date->coeff_A2 * filter_date->filter_W2;
-
-
//
-
// Y(0) = Gain * (B0 * W0 + B1 * W1 + B2 * W2)
-
//
-
filter_date->filter_out = filter_date->coeff_GAIN * ( filter_date->coeff_B0 * filter_date->filter_W0 \
-
+ filter_date->coeff_B1 * filter_date->filter_W1 \
-
+ filter_date->coeff_B2 * filter_date->filter_W2 );
-
-
//
-
// Update last data
-
//
-
filter_date->filter_W2 = filter_date->filter_W1;
-
filter_date->filter_W1 = filter_date->filter_W0;
-
-
//
-
// Return Value
-
//
-
return(filter_date->filter_out);
-
}
仿真结果：
与方法一完全一样。
参考文献：
1、Simulink 窄带陷波滤波器（Notch filter)仿真到代码生成
2、温故知新（五）——三参数陷波滤波器离散化推导及MATLAB实现
3、陷波器及其算法（基于C语言）
PLECS仿真：notch_filter.plecs
1076
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
