華北電力大學(xué)數(shù)字信號處理實(shí)驗(yàn)六

1、實(shí)驗(yàn)六 IIR數(shù)字濾波器設(shè)計(jì)及應(yīng)用一：實(shí)驗(yàn)?zāi)康募由罾斫釯IR數(shù)字濾波器的特性，掌握IIR數(shù)字濾波器的設(shè)計(jì)原理與設(shè)計(jì)方法，以及IIR數(shù)字濾波器的應(yīng)用。二：實(shí)驗(yàn)原理N階IIR數(shù)字濾波器的系統(tǒng)函數(shù)為: IIR數(shù)字濾波器的設(shè)計(jì)主要通過成熟的模擬濾波器設(shè)計(jì)方法來實(shí)現(xiàn)：將數(shù)字濾波器設(shè)計(jì)指標(biāo)轉(zhuǎn)換為模擬濾波器設(shè)計(jì)指標(biāo)，設(shè)計(jì)出相應(yīng)的模擬濾波器H(s)，再經(jīng)過脈沖響應(yīng)不變法或雙線性變換法得到所需的IIR數(shù)字濾波器H(z)。 IIR數(shù)字濾波器設(shè)計(jì)的重要環(huán)節(jié)是模擬原型低通濾波器的設(shè)計(jì)，主要包括Butterworth、Chebyshev和橢圓等濾波器。MATLAB 信號處理工具箱中提供了IIR濾波器設(shè)計(jì)的函數(shù)。IIR

2、濾波器階數(shù)選擇 buttord -巴特沃斯(Butterworth)濾波器階數(shù)選擇。 cheb1ord -切比雪夫(Chebyshev)I 型濾波器階數(shù)選擇。 cheb2ord -切比雪夫(Chebyshev)II 型濾波器階數(shù)選擇。 ellipord -橢圓(Elliptic)濾波器階數(shù)選擇。IIR 濾波器設(shè)計(jì) butter -巴特沃斯(Butterworth)濾波器設(shè)計(jì) cheby1 -切比雪夫(Chebyshev)I 型濾波器設(shè)計(jì) cheby2 -切比雪夫(Chebyshev)II 型濾波器設(shè)計(jì) ellip -橢圓(Elliptic)濾波器設(shè)計(jì) maxflat -通用的巴特沃斯(Butt

3、erworth)低通濾波器設(shè)計(jì) yulewalk -Yule-Walker 濾波器設(shè)計(jì)(直接數(shù)字濾波器設(shè)計(jì)法) 1. Butterworth濾波器設(shè)計(jì) Butterworth濾波器是通帶、阻帶都單調(diào)衰減的濾波器。(1)調(diào)用buttord函數(shù)確定巴特沃斯濾波器的階數(shù)，格式為 N,Wc = buttord(Wp,Ws,Ap,As) 輸入?yún)?shù)：Ap,As為通帶最大衰減和阻帶最小衰減，以dB為單位。 Wp,Ws為歸一化通帶截頻和阻帶截頻，0<Wp,Ws<1 。輸出參數(shù)：N為濾波器的階數(shù)；Wc為截頻，0 < Wc < 1。 (2)調(diào)用butter函數(shù)設(shè)計(jì)出巴特沃斯濾波器，格式為

4、b,a = butter(N,Wc,options) 輸入?yún)?shù)：N和Wc是buttord函數(shù)返回的參數(shù)，含義見上。 Options=low, high, bandpass, stop, 分別對應(yīng)低通、高通、帶通、帶阻，默認(rèn)情況下為低通或帶通。輸出參數(shù)：b和a為設(shè)計(jì)出的IIR數(shù)字濾波器H(s)的分子多項(xiàng)式和分母多項(xiàng)式的系數(shù)矩陣。2. Chebyshev I型濾波器設(shè)計(jì)Chebyshev I型濾波器為通帶紋波控制器：在通帶呈現(xiàn)紋波特性，在阻帶單調(diào)衰減。N,Wc = cheb1ord(Wp, Ws, Ap, As)b,a = cheby1(N,Ap,Wc,options)參數(shù)含義與butter中參數(shù)

5、一致。 2. Chebyshev II 型濾波器設(shè)計(jì) Chebyshev II 型濾波器為阻帶紋波控制器： 在阻帶呈現(xiàn)紋波特性。 N,Wc = cheb2ord(Wp, Ws, Ap, As) b,a = cheby2(N,As,Wc,options)3. 橢圓濾波器設(shè)計(jì) 橢圓濾波器在通阻帶都呈現(xiàn)紋波特性。 N,Wc = ellipord(Wp,Ws,Ap,As) b,a = ellip(N,Ap,As,Wc,options)三：實(shí)驗(yàn)內(nèi)容1（1）N,Wc=buttord(0.250,0.677,3,60)b,a=butter(N,Wc)freqz(b,a);axis(0,1,-120,0);g

6、rid ontitle('巴特沃斯低通數(shù)字濾波器')（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 ontitle('巴特沃斯高通數(shù)字濾波器')（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,

7、-120,0);grid ontitle('巴特沃斯帶通數(shù)字濾波器')N = 40Wc = 0.2499 0.6701b = 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 -

8、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 th

9、rough 81 0.0000 0 -0.0000 0 0.0000 0 -0.0000 0 0.0000a = 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.0

10、794 -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.00

11、00 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

12、.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 ontitle('巴特沃斯帶阻數(shù)字濾波器')N = 40Wc = 0.2499 0.6701b = 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.000

13、0 -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

14、-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

15、.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.0000a = 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 -

16、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.

17、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.00

18、00 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.00003（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('時(shí)域波形 '); Xm=fft(x,N)/N; f=(-(N-1)/2:(N-

19、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 ontitle('巴特沃斯低通數(shù)字濾波器')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('時(shí)域波形

20、'); 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('低通濾波后時(shí)域波形')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(

21、N,Wc,'high')freqz(b,a);axis(0,1,-120,0);grid onT0=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('時(shí)域波形'); 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);subpl

22、ot(4,1,3);stem(k,y);title('高通濾波后時(shí)域波形')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('型切比雪夫?yàn)V波器1'); g

23、rid onWp2 =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('型切比雪夫?yàn)V波器2'); grid onWp3 =830 870/4000; Ws3 =830-20 870+20/4000; Rp3 = 3; Rs3 = 40;N3,Wn3=cheb1ord(Wp3,Ws3,Rp3,Rs3)b3,a3 =cheby

24、1(N3,Rp3,Wn3); figure; freqz(b3,a3,512,8000); title('型切比雪夫?yàn)V波器3'); grid onWp4 =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('型切比雪夫?yàn)V波器4'); grid on; k=0:1:500;x=sin(2/8000)*770*p

25、i*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

26、/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('型切比雪夫?yàn)V波器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); fi

27、gure; freqz(b2,a2,512,8000); title('型切比雪夫?yàn)V波器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('型切比雪夫?yàn)V波器3'); grid on;k=0:1:500; x=sin(2/8000)*770*pi*k)+sin(2/800

28、0)*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)+

29、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/

30、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,

31、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,Wn

32、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); 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

33、/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(

34、9;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);Y1

35、3=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=filt

36、er(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);

37、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

38、(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/4

39、000; 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; n

40、1,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

41、(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),'

42、;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)

43、,'g-',k,abs(fftshift(ym03),'g-'); legend('1209HZ', '1336HZ','1477HZ'); title('第二組中"0"濾波后的頻譜');Y11=filter(B1,A1,x1); Y12=filter(B2,A2,x1);Y13=filter(B3,A3,x1);Y14=filter(B4,A4,x1);Ym11=fft(Y11,N)/N;Ym12=fft(Y12,N)/N; Ym13=fft(Y13,N)/N; Ym14=fft(Y14,N)/N