华北电力大学数字信号处理实验六.doc

实验六 IIR数字滤波器设计及应用
一：实验目的
加深理解IIR数字滤波器的特性，掌握IIR数字滤波器的设计原理与设计方法，以及IIR数字滤波器的应用。

二：实验原理
N阶IIR数字滤波器的系统函数为:

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)
输入参数：ApAs为通带最大衰减和阻带最小衰减，以dB为单位。
Wp,Ws为归一化通带截频和阻带截频，0为设计出的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);
[b1,a1] =cheby1(N1,Rp1,Wn1);
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);
[b4,a4] =cheby1(N4,Rp4,Wn4);
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','1336HZ','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;
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);
Ym01=fft(Y01,N)/N;
Ym02=fft(Y02,N)/N;
Ym03=fft(Y03,N)/N;
Ym04=fft(Y04,N)/N;
figure;
plot(k,abs(fftshift(Ym01)),k,abs(fftshift(Ym02)),'y--',k,abs(fftshift(Ym03)),'g--',k,abs(fftshift(Ym04)),'r--');
legend('697HZ', '770HZ','852HZ','941HZ'); title('第一组中"0"滤波后的频谱');
y01=filter(b1,a1,x0);
y02=filter(b2,a2,x0);
y03=filter(b3,a3,x0);
ym01=fft(y01,N)/N;
ym02=fft(y02,N)/N;
ym03=fft(y03,N)/N;
figure;
plot(k,abs(fftshift(ym01)),k,abs(fftshift(ym02)),'g--',k,abs(fftshift(ym03)),'g--');
legend('1209HZ', '1336HZ','1477HZ');
title('第二组中