电力电子技术

1、实验四数字滤波器设计一、实验原理数字滤波器的设计，就是要找到滤波器的传输函数G(z)，使其频率响应Gg逼近于给 定的频率响应指标，这个过程称为数字滤波器设计。在得到G(z)之后，就以适当的滤波器结 构的形式去实现它。根据数字滤波器的类型的不同，数字滤波器可能有不同的设计方法，比如， 对于IIR数字滤波器设计来说，需要先设计所对应的模拟滤波器，然后利用传输函数的转换， 就可得到IIR数字滤波器的传输函数。模拟低通滤波器的设计模拟滤波器的设计指标设计滤波器，首先要理解滤波器的技术指标(Specification)。就模拟滤波器来讲，其技 术指标可用图4.1来说明。Transitionband图4.

2、1归一化的模拟低通滤波器的典型幅度响应指标模拟滤波器设计涉及到的技术指标主要有四个，即：通带截止频率Qp，单位为弧度/秒；阻带截止频率Qs，单位为弧度/秒；通带最大允许衰减ap，单位为dB；阻带最小衰减as，单位为dB。其中指标ap和as与图4.1中的参数之间的关系描述如下。a = 20log1o J1 + 2dB(4.1)a = 20log1o A dB(4.2)用MATLAB设计模拟滤波器的步骤1)、Butterworth 逼近a、根据设计指标，首先确定Buterworth低通滤波器的阶N和3dB截止频率Qc(wn), MATLAB 函数为N,wn = buttord(wp,ws,Rp,R

3、s, s ):其中wp、ws、Rp和Rs分别为滤波器的技术指标通带截止频率Qp、阻带截止频率Qs、通 带最大允许衰减ap和阻带最小衰减as，Qp和Qs的单位为弧度/秒 (rad/sec)，ap和as的单位 为dB。b、根据Buterworth低通滤波器的阶N和3dB截止频率Qc，确定Buterworth低通滤波器 的系统函数Ha(s)。(4.3)QnQn(4.3)Ha (S = 3 二在 N H (S - p )ll=1其中Dn(s)为Buterworth多项式。在计算出了 Buterworth低通滤波器的3-dB截止频率和 阶N之后，系统函数可利用MATLAB函数来求得：b,a = butt

4、er(N,wn,s);其中wn为3-dB截止频率,b和a分别为系统函数的分子和分母多项式系数向量。c、根据系统函数，计算并绘制Buterworth低通滤波器的频率响应特性曲线，验证滤波器 的频率响应是否满足技术指标的要求。相关函数：H,w = freqs(b,a)计算滤波器的频率响应H，w为返回的与H相对应的频率点，单位为弧度/秒(rad/sec)。plot(w/(2*pi),20*log(abs(h)绘制滤波器的幅度增益特性曲线。2)、Chebyshev 1 型逼近a、根据设计指标，首先确定Chebyshev低通滤波器的阶N和3dB截止频率Qc(wn)，MATLAB 函数为N,wn = ch

5、eby1ord(wp,ws,Rp,Rs, s ):其中wp、ws、Rp和Rs分别为滤波器的技术指标通带截止频率Qp、阻带截止频率Qs、通 带最大允许衰减ap和阻带最小衰减as，Qp和Qs的单位为弧度/秒 (rad/sec)，ap和as的单位 为dB。b、根据Chebyshev低通滤波器的阶N和3dB截止频率Qc，确定Buterworth低通滤波器的 系统函数Ha(s)。b,a= cheby1(N,Rp,wn,s);其中wn为3-dB截止频率,b和a分别为系统函数的分子和分母多项式系数向量。c、根据系统函数，计算并绘制Chebyshev低通滤波器的频率响应特性曲线，验证滤波器的 频率响应是否满足

6、技术指标的要求。这一步骤与Butterworth低通滤波器的步骤是一样的。2、数字滤波器的设计数字滤波器的设计指标设计滤波器，首先要理解滤波器的技术指标(Specification)。数字低通滤波器的技术指 标如图4.2来说明。Pass b midStop Pass b midStop bn ndTinilfiltLOLl图4.2归一化的数字低通滤波器的典型幅度响应指标将图4.2与图4.1比较，发现不管是模拟滤波器还是数字滤波器，其技术指标是一样的。 如果有不同的话，那就是数字频率的最高值是n1)、IIR数字低通滤波器的设计IIR数字低通滤波器的设计方法主要有两种，一种是所谓的冲激响应不变法

7、(Impulse-Invariance Method)，另一种就是双线性变换法(Bilinear transfromationmethod)。 本实验只采用双线性变换的方法设计IIR数字低通滤波器。设计步骤如下第一步：根据数字滤波器的技术指标，利用双线性变换的方法，将数字滤波器的技术指标 转换成对应的模拟低通滤波器(模拟滤波器原型)的技术指标，实际上只转换相应的截止频率 即可，通带衰减和阻带衰减不变。转换公式为O。=tan(一)(4.4)2第二步:按照前述的模拟低通滤波器的设计方法，设计出模拟低通滤波器的系统函数H (s)。a第三步：利用下面的公式，将模拟低通滤波器原型的系统函数转换为所设计的

8、数字低通滤波器的系统函数：1 + sz =(4.5)1 - s2)、FIR数字低通滤波器的设计FIR数字滤波器的设计方法有多种，本实验只要求掌握窗函数法( Windowed Fourier Series)。其设计原理是，根据理想低通数字滤波器的单位冲激响应hjn，通过给加窗和移位 的方法得到一个长度为有限长的序列hn，并将这个序列看作是我们所需要设计的FIR滤波器 的单位冲激响应。由于hn是由hdn得到的，因此，hn的频率响应特性的形状一定会与hjn 的频率响应特性有相似之处。为了使我们所设计的FIR滤波器符合相应的技术指标，加窗时，可以选择不同的窗函数， 如Hanning窗、Hamming窗

9、和Blackman窗等。不同的窗函数，对FIR滤波器的指标的影响是不 同的。利用窗函数法设计FIR数字滤波器的步骤如下第一步：确定FIR数字滤波器的技术指标。第二步：确定理想低通数字滤波器的截止频率3 c，通常按照下面的公式确定:(4.6)其中，3 p和3 s分别是FIR滤波器的通带截止频率和阻带截止频率。 第三步：根据阻带衰减ap，选择合适的窗函数wn。第四步：计算FIR滤波器的过渡带宽AW=W? -W，并按照下式计算窗函数的长度，也就是FIR滤波器的长度:(4.7) c M (4.7)Aw其中c为不同窗函数所对应的过渡带宽系数，可查教材P454的表7.2得到。第五步：按照下式确定FIR滤波

10、器的单位冲激响应： / / M 一 I、(4.8)烦=岫气2| ).心(4.8)M 1 兀(n )2第六步：验证所设计的FIR滤波器是否满足技术指标的要求。按照上面的六个步骤，所设计的FIR滤波器一般是能够满足技术指标要求的。但也可能有 个别指标不能满足，比如，实际的通带截止频率与指标3 p可能不一致，这时可以适当调整 c之值，然后重复第四、五、六步。二、实验内容实验前，必须首先阅读本实验原理，了解所给的MATLAB相关函数，并根据实验内容的要求， 事先编写好实验用MATLAB程序。Project 1：模拟低通模波器设计(Analog lowpass filter design)Q4-l使用M

11、ATALB内部函数，编写MATLAB程序Q4_1，设计满足技术指标wp,ws, ap, as要求的Buterworth模拟低通滤波器。要求1、程序中，各项技术指标从键盘上键入；2、计算并绘制出滤波器的单位冲激响应h(t)，滤波器的幅度和相位频率响应曲线以及滤 波器的幅度增益Gain曲线。3、执行程序后，能够在命令窗口给出所设计的模拟低通滤波器的传递函数Ha(s)分子和分 母多项式系数向量。抄写程序 Q4_1如下：%Q4_1 clear,close all; wp = 2*pi*2000; ws = 2*pi*5000;Rp = 1; Rs = 50; N,wn = buttord(wp,ws,

12、Rp,Rs,s); b,a = butter(N,wn,s); H,w = freqs(b,a); ht=impulse(b,a); subplot(2,2,1) plot(w/(2火pi),abs(H);axis(0,6000,min(abs(H),max(abs(H) title(magnitude responszse of the butterworth filter) xlabel(Frequency index w/2火pi) subplot(2,2,2) plot(w/(2火pi),angle(H); axis(0,6000,min(angle(H),max(angle(H);t

13、itle(frequency response of the butterworth filter) xlabel(Frequency index w/2火pi) subplot(2,2,3) plot(w/(2火pi),20*log(abs(H);axis(0,6000,-150,10)title(The gain of the butterworth filter) xlabel(Frequency index w/2火pi) subplot(2,2,4) plot(ht) %axis(0,1,min(ht),max(ht) grid on disp(The numerator coeff

14、icient vector:) bdisp(The denominator coefficient vector:) aQ4-2给定技术指标如下：通带截止频率 Qp = 2nX 2000Hz rad/s；阻带截止频率 Qs = 2nX 5000Hz rad/s；通带最大允许衰减ap = 1dB；阻带最小衰减as = 50dB。执行程序 Q4_1,得到的模拟低通滤波器的单位冲激响应和其频率响应曲线图如下:magnitude responszse of the butterworth filter frequency response of the butterworth filter020004

15、0006000Frequency index w/2*pi30200040006000Frequency index w/2*pi210-1-2Frequency index w/2*piThe gain of the butterworth filterFrequency index w/2*pi-3Frequency index w/2*piThe gain of the butterworth filterFrequency index w/2*pi根据在命令窗口显示的滤波器的传递函数的分子和分母多项式系数，写出你所设计的模拟低 通滤波器的传递函数：Ha( s)从你设计的滤波起的幅度增益

16、特性曲线上，可以看出，通带截止频率为Qp =2nX2000Hz rad/s时，增益G(Qp)=阻带截止频率Qs =2nX 5000Hz rad/s时，增益G(Qs)=从而可以得出结论，你所设计的滤波器是/否满足技术指标的要求。Project 2：数字滤波器设计(Digital filter design)Q4-3 使用MATALB内部函数，编写MATLAB程序Q4_3，设计满足技术指标wp,ws, ap, as要 求的IIR数字低通滤波器，模拟滤波器原型(Analog prototype)为Buterworth模拟低通滤波器。 要求1、程序中，各项技术指标从键盘上键入；2、计算并绘制出滤波器的

17、单位冲激响应hn，滤波器的幅度和相位频率响应曲线以及滤 波器的幅度增益Gain曲线。3、执行程序后，能够在命令窗口给出所设计的IIR滤波器的传递函数H(z)分子和分母多 项式系数向量。抄写程序 Q4_3如下：%Q4_3 clear,close all; wp = 0.2*pi; ws = 0.3火pi;Wp = tan(wp/2);Ws = tan(ws/2);Rp = 1;Rs = 30;N,wn = buttord(Wp,Ws,Rp,Rs,s);b,a= butter(N,wn,s);b,a = bilinear(b,a,1/2);H,w = freqz(b,a);subplot(211)

18、plot(w/(pi),20*log(abs(H); axis(0,0.5,-150,0) grid ontitle(The gain of the IIR filter) xlabel(Frequency index w/pi) subplot(212) h = impz(b,a,50); n = 0:49;stem(n,h,.) title(The impulse response)Q4-4给定技术指标如下：通带截止频率W p = 0.2n rad；阻带截止频率W s = 0.3n rad；通带最大允许衰减ap = 1dB；阻带最小衰减as = 30dB。执行程序 Q4_3，得到的IIR滤

19、波器的单位冲激响应和其频率响应曲线图如下：The gain of the IIR filter-50;-100-150 TOC o 1-5 h z 00.050.10.150.20.250.30.350.40.450.5Frequency index w/piThe impulse response 0.3 I0.2 -JT-0.1 -. T-0心.-0.111|1111105101520253035404550根据在命令窗口显示的滤波器的传递函数的分子和分母多项式系数，写出你所设计的IIR滤波器的传递函数：H (z)从你设计的滤波起的幅度增益特性曲线上，可以看出，通带截止频率为W p = 0

20、.2n rad时，增益G(w p) = dB阻带截止频率W s = 0.3n rad时，增益G(w s) = dB从而可以得出结论，你所设计的IIR滤波器是/否满足技术指标的要求。Q45给定一个心电图的抽样信号：xn = -4, -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66,-32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0

21、, -2, -2, -2, -2, 0。编写MATLAB程序Q4_5,利用Q4-4题中设计的IIR数字滤波器，对信号xn进行滤波处理，得到滤波后的序列yn。要求使用两幅图形：在第一幅图形窗口中分别绘制滤波器的输入信号xn的时域波形和滤波器的输出信号yn的时域波形；在第二幅图形窗口中分别绘制滤波器的输入信号xn、输出信号y n和单位冲激响应hn的DFT的幅度；抄写程序Q4_5如下：%Q4_5 clear,close all;x = -4, -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -

22、60, -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0, -2, -2, -2, -2, 0; dw = 0.01;w = -pi:dw:pi; L = length(x); n = 0:L-1; subplot(222) stem(n,x,.) title(The original signal) grid %= wp = 0.2*pi; ws = 0.3火pi; Wp = tan(wp/2); Ws = tan(ws/2); Rp

23、= 1; Rs = 30; N,wn = buttord(Wp,Ws,Rp,Rs,s); b,a= butter(N,wn,s); b,a = bilinear(b,a,0.5) H = freqz(b,a,w); %= subplot(221) plot(w/pi,20*log(abs(H); axis(0,0.5,-150,0) grid on title(The gain of the IIR filter) xlabel(Frequency index w/pi) %= subplot(223) h = impz(b,a,L); stem(n,h, .) title(The impul

24、se response) grid xlabel(Time index n) %= X = freqz(x,1,w);Y = X.*H;y = dw*Y*exp(j*w*n)/(2*pi);%= subplot(224)stem(n,y,.*)title(* The output signal of the filter) gridxlabel(Time index n*) %= figure1 = 0 : 2*L-2;subplot(224) yl=conv(x,h) stem(lr ylr * .) axis ( 0AL-l,-100r50) %=subplot(221)X=fft(x,L

25、)stem(n,X);title(XK)%= subplot(222)H=fft(h,L)stem(n,H);title(*HK)%= subplot(223)Y=fft(y,L)stem(n,Y);title(YK1)执行 程序Q4_5,得到信号的各种图形如下The original signalThe output signal of the filterXKHK204060YKQ4-6使用窗函数法，按照前面给出的窗函数法设计步骤，编写MATLAB程序Q4_6，设计满 足技术指标wp,ws, ap, as要求的FIR数字低通滤波器，要求1、程序中，各项技术指标从键盘上键入；2、计算并绘制出

26、滤波器的单位冲激响应hn，滤波器的幅度和相位频率响应曲线以及滤 波器的幅度增益Gain曲线。3、执行程序后，能够在命令窗口给出所设计的FIR滤波器的传递函数H(z)或单位冲激响 应 hn。抄写程序 Q4_6如下：% Q4_6 clear,close all; wp=input(Type in the passband edge frequency (in Radian) wp=); ws=input(Type in the passband edge frequency (in Radian) ws=); ap = input(Type in the passband attenuation

27、ap =:); as = input(Type in the passband attenuation as =:); deltaw = ws-wp;wc=(wp+ws)/2;%= if as 75.3M = ceil(5.56*pi/deltaw)L = 2*M+1;n=0:L-1;win=blackman(L);elseif as 54.5M = ceil(3.32*pi/deltaw)L = 2*M+1;n=0:L-1;win=hamming(L);elseif as 43.9M = ceil(3.11*pi/deltaw)L = 2*M+1;n=0:L-1;win=hanning(L)

28、;elseM = ceil(0.92*pi/deltaw)L = 2*M+1;n=0:L-1;win=u(n)-u(n-L);end%= h=(wc/pi)*sinc(wc*(n-M)/pi).*win; % Impulse response %=N=2000;H=fft(h,N);H=H/max(abs(H);W=0:2*pi/N:2*pi-2*pi/N;%=subplot(222)plot(W/pi,abs(H),title(The magnitude response of the FIR filter) axis(0,1,0,max(abs(H),grid on; %= subplot

29、(223)plot(W/pi,angle(H)axis(0,1,-pi,pi),grid on;title(The phase response of the FIR filter),xlabel(Frequency (*pi) %= subplot(224)plot(W/pi,20*log10(abs(H);axis(0,1,-100,0),grid on;title(The frequency response in dB), xlabel(Frequency (*pi) %= subplot(221) stem(n,h,.), axis(0,L-1,min(h),max(h),grid

30、on;title(The impulse response of the FIR filter) xlabel(Index n) %=Q4-7给定技术指标如下：通带截止频率W p = 0.2n rad；阻带截止频率W s = 0.3n rad；通带最大允许衰减ap = 1dB；阻带最小衰减as = 30dB。执行程序 Q4_6，得到的FIR滤波器的单位冲激响应和其频率响应曲线图如下:The impulse response of the FIR filterThe magnitude response of the FIR filterThe phase response of the FIR filterThe frequency response in dBThe impulse response of the FIR filterThe magnitude response of the FIR filterThe phase response of the FIR filterThe frequency response in dB根据在命令窗口显示的滤波器的单位冲激响应，写出你所设计的FIR滤波器的传递函数:从你设计的滤波起的幅度增益特性曲线上，可以看出，通带截止频率为W p = 0.2n rad时，