陷波器设计

陷波器设计

陷波器是无限冲击响应(IIR)数字滤波器，该滤波器可以用以下常系数线性差分方程表示：

∑∑==---=M i N

i i i i n y b i n x a n y 01)()()( (1)

式中： x(n)和y(n)分别为输人和输出信号序列；i a 和i b 为滤波器系数。 对式(1)两边进行z 变换，得到数字滤波器的传递函数为：

∏∏∑∑===-=---==N i i M i i N i i i

M i i

i p z z z z b z a z H 1100)()()( (2)

式中：i z 和i p 分别为传递函数的零点和极点。

由传递函数的零点和极点可以大致绘出频率响应图。在零点处，频率响应出现极小值；在极点处，频率响应出现极大值。因此可以根据所需频率响应配置零点和极点，然后反向设计带陷数字滤波器。考虑一种特殊情况，若零点i z 在第1象限单位圆上，极点i p 在单位圆内靠近零点的径向上。为了防止滤波器系数出现复数，必须在z 平面第4象限对称位置配置相应的共轭零点*i z 、共轭极点*i p 。

这样零点、极点配置的滤波器称为单一频率陷波器，在频率ωo 处出现凹陷。而把极点设置在零的的径向上距圆点的距离为l-μ处，陷波器的传递函数为：

)

)1()()1(())(()(2121z z z z z z z z z H μμ------= (3) 式(3)中μ越小，极点越靠近单位圆，则频率响应曲线凹陷越深，凹陷的宽度也越窄。当需要消除窄带干扰而不能对其他频率有衰减时，陷波器是一种去除窄带干扰的理想数字滤波器。

当要对几个频率同时进行带陷滤波时，可以按(2)式把几个单独频率的带陷滤波器(3)式串接在一起。

一个例子：设有一个输入，它由50Hz 信号和100Hz 信号组成。50Hz 是一个干扰信号，要设计一个50 Hz 的带陷滤波器，采样频率为400Hz 。

4/400/5021ππω=⨯=

因此z 平面上的零极点可设置为

4/14

/1999.0ππj j e

p e z ±±== 展开式为

7063

7064)707.0707.0(999.0)4sin 4(cos

999.0999.0707.0707.02

2224sin 4cos 4/14/1j j j e p j j j e z j j ±=±=±=±=±=±=±±π

ππ

π

ππ＝＝ 它的传递函数为

2

12

1221111999.04126.11414.11999.04126.11414.1)7063.07063.0)(7063.07063.0()

707.0707.0)(707.0707.0()

)(())(()(----**+-+-=+-+-=+---+---=----=z z z z z z z z j z j z j z j z p z p z z z z z z H

因此分子系数是[1 1.414 1]；分母系数是[1 1.4126 0.999]。

差分方程有

)

2()3()1()2()2()3()1()2()()()

2()3()1()2()()2()3()1()2()(-----+-+=∴-+-+=-+-+n y a n y a n x b n x b n x n y n x b n x b n x n y a n y a n y

程序清单有

% 50Hz notch filter

% sample frequency=400

%

clear all;

clc;

b=[1 -sqrt(2) 1];

a=[1 -sqrt(2)*0.999 0.999];

[db, mag, pha, grd, w]=freqz_m(b, a); subplot(221); plot(w*200/pi, db); title(' Magnitude Response' ); xlabel('frequency in Hz'); ylabel('dB'); axis([0, 100, -200, 5]); set(gca, 'XTickMode', 'manual', 'XTick', [0, 50, 100]);

set(gca, 'YTickmode', 'manual', 'YTick', [-200, -100, 0]); grid title('Notch filter response');

t0=1:8000;

t=1:256;

t1=1:100;

t2=1:128;

x=sin(2*pi*50*t0/400)+0.5*sin(2*pi*100*t0/400);

x1=x(t);

y=filter(b,a,x1);

subplot(222); plot(x1);

title('Original waveform');

X=fft(x1);

subplot(223); plot(t2*400/256,abs(X(t2)));

xlabel('frequency in Hz'); ylabel('|H|'); axis([0, 200, 0, 150]); title('Spectrum for original');

set(gca, 'XTickMode', 'manual', 'XTick', [0, 50, 100, 150]); set(gca, 'YTickmode', 'manual', 'YTick', [50, 100]); grid

y=filter(b,a,x);

x1=y(t+7600);

X=fft(x1);