华北电力大学数字信号处理实验六
华北电力大学数字信号处理实验六
一：实验目的
加深明白得IIR 数字滤波器的特性，把握IIR 数字滤波器的设计原理与设计方法，以及IIR 数字滤波器的应用。
二：实验原理
N 阶IIR 数字滤波器的系统函数为: i i N i j j M
j z a z b z H -=-=∑∑+=10
1)( N
N N N M M
M M z a z a z a z b z b z b b ----------++++++++=)1(111)1(11101ΛΛ
IIR 数字滤波器的设计要紧通过成熟的模拟滤波器设计方法来实现：将数字滤波器设计指标转换为模拟滤波器设计指标，设计出相应的模拟滤波器H(s)，再通过脉冲响应不变法或双线性变换法得到所需的IIR 数字滤波器H(z)。
IIR 数字滤波器设计的重要环节是模拟原型低通滤波器的设计，要紧包括Butterworth 、Chebyshev 和椭圆等滤波器。
MATLAB 信号处理工具箱中提供了IIR 滤波器设计的函数。 IIR 滤波器阶数选择
buttord -巴特沃斯(Butterworth)滤波器阶数选择。 cheb1ord -切比雪夫(Chebyshev)I 型滤波器阶数选择。 cheb2ord -切比雪夫(Chebyshev)II 型滤波器阶数选择。 ellipord -椭圆(Elliptic)滤波器阶数选择。 IIR 滤波器设计
butter -巴特沃斯(Butterworth)滤波器设计 cheby1 -切比雪夫(Chebyshev)I 型滤波器设计 cheby2 -切比雪夫(Chebyshev)II 型滤波器设计 ellip -椭圆(Elliptic)滤波器设计
maxflat -通用的巴特沃斯(Butterworth)低通滤波器设计
yulewalk -Yule-Walker 滤波器设计(直截了当数字滤波器设计法) 1. Butterworth 滤波器设计
Butterworth 滤波器是通带、阻带都单调衰减的滤波器。
(1)调用buttord函数确定巴特沃斯滤波器的阶数，格式为[N,Wc] = buttord(Wp,Ws,Ap,As)
输入参数：Ap,As为通带最大衰减和阻带最小衰减，以dB为单位。
Wp,Ws为归一化通带截频和阻带截频，0 输出参数：N为滤波器的阶数；Wc为截频，0 < Wc < 1。 (2)调用butter函数设计出巴特沃斯滤波器，格式为[b,a] = butter(N, Wc,options) 输入参数：N和Wc是buttord函数返回的参数，含义见上。 Options=’low’, ’high’, ’bandpass’, ’stop’,分不对应低通、高通、带通、带阻，默认情形下为低通或带通。 输出参数：b和a为设计出的IIR数字滤波器H(s)的分子多项式和分母多项式的系数矩阵。 2. Chebyshev I型滤波器设计 Chebyshev I型滤波器为通带纹波操纵器：在通带出现纹波特性，在阻带单调衰减。 [N,Wc] = cheb1ord(Wp, Ws, Ap, As) [b,a] = cheby1(N,Ap,Wc,options) 参数含义与butter中参数一致。 2. Chebyshev II 型滤波器设计 Chebyshev II 型滤波器为阻带纹波操纵器： 在阻带出现纹波特性。 [N,Wc] = cheb2ord(Wp, Ws, Ap, As) [b,a] = cheby2(N,As,Wc,options) 3. 椭圆滤波器设计 椭圆滤波器在通阻带都出现纹波特性。 [N,Wc] = ellipord(Wp,Ws,Ap,As) [b,a] = ellip(N,Ap,As,Wc,options) 三：实验内容 1（1） [N,Wc]=buttord(0.250,0.677,3,60) [b,a]=butter(N,Wc) freqz(b,a); axis([0,1,-120,0]); grid on title('巴特沃斯低通数字滤波器') （2） [N,Wc]=buttord(0.250,0.677,3,60) [b,a]=butter(N,Wc,'high') freqz(b,a); axis([0,1,-120,0]); grid on title('巴特沃斯高通数字滤波器') （3） Wp =[0.25 0.67]; Ws =[0.25-0.03 0.67+0.03]; Rp = 3; Rs = 60; [N,Wc]=buttord(Wp,Ws,Rp,Rs) [b,a]=butter(N,Wc,'bandpass') freqz(b,a); axis([0,1,-120,0]); grid on title('巴特沃斯带通数字滤波器') N = 40 Wc = 0.2499 0.6701 b = Columns 1 through 9 0.0000 0 -0.0000 0 0.0000 0 -0.0000 0 0.0000 Columns 10 through 18 0 -0.0000 0 0.0000 0 -0.0000 0 0.0000 0 Columns 19 through 27 -0.0001 0 0.0003 0 -0.0007 0 0.0017 0 -0.0037 Columns 28 through 36 0 0.0072 0 -0.0125 0 0.0195 0 -0.0276 0 Columns 37 through 45 0.0353 0 -0.0408 0 0.0429 0 -0.0408 0 0.0353 Columns 46 through 54 0 -0.0276 0 0.0195 0 -0.0125 0 0.0072 0 Columns 55 through 63 -0.0037 0 0.0017 0 -0.0007 0 0.0003 0 -0.0001 Columns 64 through 72 0 0.0000 0 -0.0000 0 0.0000 0 -0.0000 0 Columns 73 through 81 0.0000 0 -0.0000 0 0.0000 0 -0.0000 0 0.0000 a = 1.0e+005 * Columns 1 through 9 0.0000 -0.0001 0.0003 -0.0011 0.0030 -0.0074 0.0160 -0.0318 0.0585 Columns 10 through 18 -0.1008 0.1637 -0.2519 0.3692 -0.5174 0.6952 -0.8980 1.1176 -1.3427 Columns 19 through 27 1.5597 -1.7542 1.9125 - 2.0234 2.0794 -2.0773 2.0188 -1.9098 1.7596 Columns 28 through 36 -1.5798 1.3828 -1.1803 0.9828 -0.7985 0.6332 -0.4902 0.3705 -0.2734 Columns 37 through 45 0.1970 -0.1386 0.0953 -0.0639 0.0419 -0.0268 0.0167 -0.0102 0.0061 Columns 46 through 54 -0.0035 0.0020 -0.0011 0.0006 -0.0003 0.0002 -0.0001 0.0000 -0.0000 Columns 55 through 63 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 Columns 64 through 72 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 Columns 73 through 81 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 >> （4） Wp =[0.25 0.67]; Ws =[0.25-0.03 0.67+0.03]; Rp = 3; Rs = 60; [N,Wc]=buttord(Wp,Ws,Rp,Rs) [b,a]=butter(N,Wc,'stop') freqz(b,a); axis([0,1,-120,0]); grid on title('巴特沃斯带阻数字滤波器') N = 40 Wc = 0.2499 0.6701 b = 1.0e+005 * Columns 1 through 7 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 Columns 8 through 14 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 Columns 15 through 21 0.0001 -0.0001 0.0003 -0.0007 0.0015 -0.0029 0.0056 Columns 22 through 28 -0.0102 0.0179 -0.0305 0.0502 -0.0798 0.1227 -0.1828 Columns 29 through 35 0.2636 -0.3686 0.4998 -0.6576 0.8396 -1.0409 1.2534 Columns 36 through 42 -1.4660 1.6661 -1.8401 1.9753 -2.0610 2.0904 -2.0610 Columns 43 through 49 1.9753 -1.8401 1.6661 -1.4660 1.2534 -1.0409 0.8396 Columns 50 through 56 -0.6576 0.4998 -0.3686 0.2636 -0.1828 0.1227 -0.0798 Columns 57 through 63 0.0502 -0.0305 0.0179 -0.0102 0.0056 -0.0029 0.0015 Columns 64 through 70 -0.0007 0.0003 -0.0001 0.0001 -0.0000 0.0000 -0.0000 Columns 71 through 77 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 Columns 78 through 81 -0.0000 0.0000 -0.0000 0.0000 a = 1.0e+005 * Columns 1 through 7 0.0000 -0.0001 0.0003 -0.0011 0.0030 -0.0074 0.0160 Columns 8 through 14 -0.0318 0.0585 -0.1008 0.1637 -0.2519 0.3692 -0.5174 Columns 15 through 21 0.6952 -0.8980 1.1176 -1.3427 1.5597 -1.7542 1.9125 Columns 22 through 28 -2.0234 2.0794 -2.0773 2.0188 -1.9098 1.7596 -1.5798 Columns 29 through 35 1.3828 -1.1803 0.9828 -0.7985 0.6332 -0.4902 0.3705 Columns 36 through 42 -0.2734 0.1970 -0.1386 0.0953 -0.0639 0.0419 -0.0268 Columns 43 through 49 0.0167 -0.0102 0.0061 -0.0035 0.0020 -0.0011 0.0006 Columns 50 through 56 -0.0003 0.0002 -0.0001 0.0000 -0.0000 0.0000 -0.0000 Columns 57 through 63 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 Columns 64 through 70 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 Columns 71 through 77 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 Columns 78 through 81 -0.0000 0.0000 -0.0000 0.0000 3（1） T0=204; N=205; T=1; k=0:T0; x=sin((2/8000)*770*pi*k)+sin((2/8000)*1209*pi*k); subplot(2,1,1); stem(k,x); title('时域波形'); Xm=fft(x,N)/N; f=(-(N-1)/2:(N-1)/2)/N/T; subplot(2,1,2); stem(f,abs(fftshift(Xm))); title('频谱图'); （2） [N,Wc]=buttord(0.1925,0.30225,3,60) [b,a]=butter(N,Wc) freqz(b,a); axis([0,1,-120,0]); grid on title('巴特沃斯低通数字滤波器') T0=204; N=205; T=1; k=0:T0; x=sin((2/8000)*770*pi*k)+sin((2/8000)*1209*pi*k); subplot(4,1,1); stem(k,x); title('时域波形'); Xm=fft(x,N)/N; f=(-(N-1)/2:(N-1)/2)/N/T; subplot(4,1,2); stem(f,abs(fftshift(Xm))); title('频谱图'); y=filter(b,a,x); subplot(4,1,3); stem(k,y); title('低通滤波后时域波形') ym=fft(y,N)/N; subplot(4,1,4); stem(f,abs(fftshift(ym))); title('低通滤波后频谱图') [N,Wc]=buttord(0.1925,0.30225,3,60) [b,a]=butter(N,Wc,'high') freqz(b,a); axis([0,1,-120,0]); grid on T0=204; N=205; T=1; k=0:T0; x=sin((2/8000)*770*pi*k)+sin((2/8000)*1209*pi*k); subplot(4,1,1); stem(k,x); title('时域波形'); Xm=fft(x,N)/N; f=(-(N-1)/2:(N-1)/2)/N/T; subplot(4,1,2); stem(f,abs(fftshift(Xm))); title('频谱图'); y=filter(b,a,x); subplot(4,1,3); stem(k,y); title('高通滤波后时域波形') ym=fft(y,N)/N; subplot(4,1,4); stem(f,abs(fftshift(ym))); title('高通滤波后频谱图') （3） Wp1 =[680 720]/4000; Ws1=[650-20 720+20]/4000; Rp1 = 3; Rs1 = 40; [N1,Wn1]=cheb1ord(Wp1,Ws1,Rp1,Rs1); freqz(b1,a1,512,8000); title('Ⅰ型切比雪夫滤波器1'); grid on Wp2 =[750 790]/4000; Ws2 =[750-20 790+20]/4000; Rp2 = 3; Rs2 = 40; [N2,Wn2]=cheb1ord(Wp2,Ws2,Rp2,Rs2) [b2,a2] =cheby1(N2,Rp2,Wn2); figure; freqz(b2,a2,512,8000); title('Ⅰ型切比雪夫滤波器2'); grid on Wp3 =[830 870]/4000; Ws3 =[830-20 870+20]/4000; Rp3 = 3; Rs3 = 40; [N3,Wn3]=cheb1ord(Wp3,Ws3,Rp3,Rs3) [b3,a3] =cheby1(N3,Rp3,Wn3); figure; freqz(b3,a3,512,8000); title('Ⅰ型切比雪夫滤波器3'); grid on Wp4 =[920 960]/4000; Ws4 =[920-20 960+20]/4000; Rp4 = 3; Rs4 = 40; [N4,Wn4]=cheb1ord(Wp4,Ws4,Rp4,Rs4); figure; freqz(b4,a4,512,8000); title('Ⅰ型切比雪夫滤波器4'); grid on; k=0:1:500; x=sin((2/8000)*770*pi*k)+sin((2/8000)*1209*pi*k); y1=filter(b1,a1,x); y2=filter(b2,a2,x); y3=filter(b3,a3,x); y4=filter(b4,a4,x); figure; plot(k,y1,k,y2,'g--',k,y3,'r--',k,y4,'y--'); title('滤波后4条输出曲线') ; legend('697HZ', '770HZ','852HZ','941HZ'); （4） Wp1 =[1180 1220]/4000; Ws1 =[1180-30 1220+30]/4000; Rp1 = 3; Rs1 = 40; [N1,Wn1]=cheb1ord(Wp1,Ws1,Rp1,Rs1) [b1,a1] =cheby1(N1,Rp1,Wn1); freqz(b1,a1,512,8000); title('Ⅰ型切比雪夫滤波器1'); grid on; Wp2 =[1310 1350]/4000; Ws2 =[1310-30 1350+30]/4000; Rp2 = 3; Rs2 = 40; [N2,Wn2]=cheb1ord(Wp2,Ws2,Rp2,Rs2); [b2,a2] =cheby1(N2,Rp2,Wn2); figure; freqz(b2,a2,512,8000); title('Ⅰ型切比雪夫滤波器2'); grid on; Wp3 =[1460 1500]/4000; Ws3=[1460-30 1500+30]/4000; Rp3 = 3; Rs3 = 40; [N3,Wn3]=cheb1ord(Wp3,Ws3,Rp3,Rs3) [b3,a3] =cheby1(N3,Rp3,Wn3); figure; freqz(b3,a3,512,8000); title('Ⅰ型切比雪夫滤波器3'); grid on; k=0:1:500; x=sin((2/8000)*770*pi*k)+sin((2/8000)*1209*pi*k); y1=filter(b1,a1,x); y2=filter(b2,a2,x); y3=filter(b3,a3,x); figure; plot(k,y1,k,y2,'g--',k,y3,'y--'); title('输出曲线') ; legend('1209HZ','1336 HZ','1477HZ'); (5) k=0:1:500; x0=sin((2/8000)*941*pi*k)+sin((2/8000)*1336*pi*k); x1=sin((2/8000)*697*pi*k)+sin((2/8000)*1209*pi*k); x2=sin((2/8000)*697*pi*k)+sin((2/8000)*1336*pi*k); x3=sin((2/8000)*697*pi*k)+sin((2/8000)*1477*pi*k); x4=sin((2/8000)*770*pi*k)+sin((2/8000)*1209*pi*k); x5=sin((2/8000)*770*pi*k)+sin((2/8000)*1336*pi*k); x6=sin((2/8000)*770*pi*k)+sin((2/8000)*1477*pi*k); x7=sin((2/8000)*852*pi*k)+sin((2/8000)*1209*pi*k); x8=sin((2/8000)*852*pi*k)+sin((2/8000)*1336*pi*k); x9=sin((2/8000)*852*pi*k)+sin((2/8000)*1477*pi*k); Wp1 =[680 720]/4000; Ws1=[650-20 720+20]/4000; Rp1 = 3; Rs1 = 40; [N1,Wn1]=cheb1ord(Wp1,Ws1,Rp1,Rs1); [B1,A1] =cheby1(N1,Rp1,Wn1); Wp2 =[750 790]/4000; Ws2 =[750-20 790+20]/4000; Rp2 = 3; Rs2 = 40; [N2,Wn2]=cheb1ord(Wp2,Ws2,Rp2,Rs2); [B2,A2] =cheby1(N2,Rp2,Wn2); Wp3 =[830 870]/4000; Ws3 =[830-20 870+20]/4000; Rp3 = 3; Rs3 = 40; [N3,Wn3]=cheb1ord(Wp3,Ws3,Rp3,Rs3); [B3,A3] =cheby1(N3,Rp3,Wn3); Wp4 =[920 960]/4000; Ws4 =[920-20 960+20]/4000; Rp4 = 3; Rs4 = 40; [N4,Wn4]=cheb1ord(Wp4,Ws4,Rp4,Rs4); [B4,A4] =cheby1(N4,Rp4,Wn4); wp1 =[1180 1220]/4000; ws1 =[1180-30 1220+30]/4000; rp1 = 3; rs1 = 40; [n1,wn1]=cheb1ord(wp1,ws1,rp1,rs1); [b1,a1] =cheby1(n1,rp1,wn1); wp2 =[1310 1350]/4000; ws2 =[1310-30 1350+30]/4000; rp2 = 3; rs2 = 40; [n2,wn2]=cheb1ord(wp2,ws2,rp2,rs2); [b2,a2] =cheby1(n2,rp2,wn2); wp3 =[1460 1500]/4000; ws3=[1460-30 1500+30]/4000; rp3 = 3; rs3 = 40; [n3,wn3]=cheb1ord(wp3,ws3,rp3,rs3); [b3,a3]=cheby1(n3,rp3,wn3); Y01=filter(B1,A1,x0); Y02=filter(B2,A2,x0); Y03=filter(B3,A3,x0); Y04=filter(B4,A4,x0); figure ; subplot(2,1,1); plot(k,Y01,k,Y02,'y--',k,Y03,'r--',k,Y04,'g--'); title('输出曲线1'); legend('697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y01=filter(b1,a1,x0); y02=filter(b2,a2,x0); y03=filter(b3,a3,x0); plot(k,y01,k,y02,'g--',k,y03,'r--'); legend('1209HZ', '1336HZ','1477HZ'); Y11=filter(B1,A1,x1); Y12=filter(B2,A2,x1); Y13=filter(B3,A3,x1); Y14=filter(B4,A4,x1); figure;subplot(2,1,1); plot(k,Y11,k,Y12,'y--',k,Y13,'g--',k,Y14,'r--'); title('输出曲线2') ; legend('697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y11=filter(b1,a1,x1); y12=filter(b2,a2,x1); y13=filter(b3,a3,x1); plot(k,y11,k,y12,'y--',k,y13,'g--'); legend('1209HZ', '1336HZ','1477HZ'); (5) T0=500; N=501; T=1;k=0:T:T0;x0=sin((2/8000)*941*pi*k)+sin((2/8000)*1336*pi*k); x1=sin((2/8000)*697*pi*k)+sin((2/8000)*1209*pi*k); x2=sin((2/8000)*697*pi*k)+sin((2/8000)*1336*pi*k); x3=sin((2/8000)*697*pi*k)+sin((2/8000)*1477*pi*k); x4=sin((2/8000)*770*pi*k)+sin((2/8000)*1209*pi*k); x5=sin((2/8000)*770*pi*k)+sin((2/8000)*1336*pi*k); x6=sin((2/8000)*770*pi*k)+sin((2/8000)*1477*pi*k); x7=sin((2/8000)*852*pi*k)+sin((2/8000)*1209*pi*k); x8=sin((2/8000)*852*pi*k)+sin((2/8000)*1336*pi*k); x9=sin((2/8000)*852*pi*k)+sin((2/8000)*1477*pi*k); Wp1 =[680 720]/4000; Ws1=[650-20 720+20]/4000; Rp1 = 3; Rs1 = 40; [N1,Wn1]=cheb1ord(Wp1,Ws1,Rp1,Rs1); [B1,A1] =cheby1(N1,Rp1,Wn1); Wp2 =[750 790]/4000; Ws2 =[750-20 790+20]/4000; Rp2 = 3; Rs2 = 40; [N2,Wn2]=cheb1ord(Wp2,Ws2,Rp2,Rs2); [B2,A2] =cheby1(N2,Rp2,Wn2); Wp3 =[830 870]/4000; Ws3 =[830-20 870+20]/4000; Rp3 = 3;