精编华北电力大学数字信号处理实验六.docx
- 文档编号:5369423
- 上传时间:2022-12-15
- 格式:DOCX
- 页数:19
- 大小:18.82KB
精编华北电力大学数字信号处理实验六.docx
《精编华北电力大学数字信号处理实验六.docx》由会员分享,可在线阅读,更多相关《精编华北电力大学数字信号处理实验六.docx(19页珍藏版)》请在冰豆网上搜索。
精编华北电力大学数字信号处理实验六
(精编)华北电力大学数字信号处理实验六
实验六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)
输入参数:
Ap,As为通带最大衰减和阻带最小衰减,以dB为单位。
Wp,Ws为归一化通带截频和阻带截频,0 输出参数: N为滤波器的阶数;Wc为截频,0 (2)调用butter函数设计出巴特沃斯滤波器,格式为[b,a]=butter(N,Wc,options) 输入参数: N和Wc是buttord函数返回的参数,含义见上。 Options=’low’,’high’,’bandpass’,’stop’,分别对应低通、高通、带通、带阻,默认情况下为低通或带通。 输出参数: b和a为设计出的IIR数字滤波器H(s)的分子多项式和分母多项式的系数矩阵。 2.ChebyshevI型滤波器设计 ChebyshevI型滤波器为通带纹波控制器: 在通带呈现纹波特性,在阻带单调衰减。 [N,Wc]=cheb1ord(Wp,Ws,Ap,As) [b,a]=cheby1(N,Ap,Wc,options) 参数含义与butter中参数一致。 2.ChebyshevII型滤波器设计 ChebyshevII型滤波器为阻带纹波控制器: 在阻带呈现纹波特性。 [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]); gridon 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]); gridon title('巴特沃斯高通数字滤波器') (3) Wp=[0.250.67];Ws=[0.25-0.030.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]); gridon title('巴特沃斯带通数字滤波器') N= 40 Wc= 0.24990.6701 b= Columns1through9 0.00000-0.000000.00000-0.000000.0000 Columns10through18 0-0.000000.00000-0.000000.00000 Columns19through27 -0.000100.00030-0.000700.00170-0.0037 Columns28through36 00.00720-0.012500.01950-0.02760 Columns37through45 0.03530-0.040800.04290-0.040800.0353 Columns46through54 0-0.027600.01950-0.012500.00720 Columns55through63 -0.003700.00170-0.000700.00030-0.0001 Columns64through72 00.00000-0.000000.00000-0.00000 Columns73through81 0.00000-0.000000.00000-0.000000.0000 a= 1.0e+005* Columns1through9 0.0000-0.00010.0003-0.00110.0030-0.00740.0160-0.03180.0585 Columns10through18 -0.10080.1637-0.25190.3692-0.51740.6952-0.89801.1176-1.3427 Columns19through27 1.5597-1.75421.9125-2.02342.0794-2.07732.0188-1.90981.7596 Columns28through36 -1.57981.3828-1.18030.9828-0.79850.6332-0.49020.3705-0.2734 Columns37through45 0.1970-0.13860.0953-0.06390.0419-0.02680.0167-0.01020.0061 Columns46through54 -0.00350.0020-0.00110.0006-0.00030.0002-0.00010.0000-0.0000 Columns55through63 0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000 Columns64through72 -0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.0000 Columns73through81 0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000 >> (4) Wp=[0.250.67]; Ws=[0.25-0.030.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]); gridon title('巴特沃斯带阻数字滤波器') N= 40 Wc= 0.24990.6701 b= 1.0e+005* Columns1through7 0.0000-0.00000.0000-0.00000.0000-0.00000.0000 Columns8through14 -0.00000.0000-0.00000.0000-0.00000.0000-0.0000 Columns15through21 0.0001-0.00010.0003-0.00070.0015-0.00290.0056 Columns22through28 -0.01020.0179-0.03050.0502-0.07980.1227-0.1828 Columns29through35 0.2636-0.36860.4998-0.65760.8396-1.04091.2534 Columns36through42 -1.46601.6661-1.84011.9753-2.06102.0904-2.0610 Columns43through49 1.9753-1.84011.6661-1.46601.2534-1.04090.8396 Columns50through56 -0.65760.4998-0.36860.2636-0.18280.1227-0.0798 Columns57through63 0.0502-0.03050.0179-0.01020.0056-0.00290.0015 Columns64through70 -0.00070.0003-0.00010.0001-0.00000.0000-0.0000 Columns71through77 0.0000-0.00000.0000-0.00000.0000-0.00000.0000 Columns78through81 -0.00000.0000-0.00000.0000 a= 1.0e+005* Columns1through7 0.0000-0.00010.0003-0.00110.0030-0.00740.0160 Columns8through14 -0.03180.0585-0.10080.1637-0.25190.3692-0.5174 Columns15through21 0.6952-0.89801.1176-1.34271.5597-1.75421.9125 Columns22through28 -2.02342.0794-2.07732.0188-1.90981.7596-1.5798 Columns29through35 1.3828-1.18030.9828-0.79850.6332-0.49020.3705 Columns36through42 -0.27340.1970-0.13860.0953-0.06390.0419-0.0268 Columns43through49 0.0167-0.01020.0061-0.00350.0020-0.00110.0006 Columns50through56 -0.00030.0002-0.00010.0000-0.00000.0000-0.0000 Columns57through63 0.0000-0.00000.0000-0.00000.0000-0.00000.0000 Columns64through70 -0.00000.0000-0.00000.0000-0.00000.0000-0.0000 Columns71through77 0.0000-0.00000.0000-0.00000.0000-0.00000.0000 Columns78through81 -0.00000.0000-0.00000.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]); gridon 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]); gridon 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=[680720]/4000; Ws1=[650-20720+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'); gridon Wp2=[750790]/4000; Ws2=[750-20790+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'); gridon Wp3=[830870]/4000; Ws3=[830-20870+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'); gridon Wp4=[920960]/4000; Ws4=[920-20960+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'); gridon; 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=[11801220]/4000; Ws1=[1180-301220+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'); gridon; Wp2=[13101350]/4000; Ws2=[1310-301350+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'); gridon; Wp3=[14601500]/4000; Ws3=[1460-301500+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'); gridon; 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=[680720]/4000; Ws1=[650-20720+20]/4000; Rp1=3; Rs1=40; [N1,Wn1]=cheb1ord(Wp1,Ws1,Rp1,Rs1); [B1,A1]=cheby1(N1,Rp1,Wn1); Wp2=[750790]/4000; Ws2=[750-20790+20]/4000; Rp2=3; Rs2=40; [N2,Wn2]=cheb1ord(Wp2,Ws2,Rp2,Rs2); [B2,A2]=cheby1(N2,Rp2,Wn2); Wp3=[830870]/4000; Ws3=[830-20870+20]/4000; Rp3=3; Rs3=40; [N3,Wn3]=cheb1ord(Wp3,Ws3,Rp3,Rs3); [B3,A3]=cheby1(N3,Rp3,Wn3); Wp4=[920960]/4000; Ws4=[920-20960+20]/4000; Rp4=3; Rs4=40; [N4,Wn4]=cheb1ord(Wp4,Ws4,Rp4,Rs4); [B4,A4]=cheby1(N4,Rp4,Wn4); wp1=[11801220]/4000; ws1=[1180-301220+30]/4000; rp1=3; rs1=40; [n1,wn1]=cheb1ord(wp1,ws1,rp1,rs1); [b1,a1]=cheby1(n1,rp1,wn1); wp2=[13101350]/4000; ws2=[1310-301350+30]/4000; rp2=3; rs2=40; [n2,wn2]=cheb1ord(wp2,ws2,rp2,rs2); [b2,a2]=cheby1(n2,rp2,wn2); wp3=[14601500]/4000; ws3=[1460-301500+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;sub
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 精编 华北电力 大学 数字信号 处理 实验