数字信号处理实验.docx

1、实验七 一、实验目的 加深理解 IIR 数字滤波器的特性，掌握 IIR 数字滤波器的设计原理与设计方法，以及 IIR数字滤波器的应用。 二、 实验原理 N 阶 IIR 数字滤波器的系统函数为 : IIR 数字滤波器的设计主要通过成熟的模拟滤波器设计方法来实现：将数字滤波器设计指标转换为模拟滤波器设计指标，设计出相应的模拟滤波器 H(s)，再经过脉冲响应不变法或双线性变换法得到所需的 IIR 数字滤波器 H(z)。 IIR 数字滤波器设计的重要环节是模拟原型低通滤波器的设计，主要包括Butterworth、 Chebyshev 和椭圆等滤波器 MATLAB 信号处理工具箱中提供了 IIR 滤波器

2、设计的函数。 IIR 滤波器阶数选择 buttord -巴特沃斯 (Butterworth)滤波器阶数选择。 cheb1ord -切比雪夫 (Chebyshev)I 型滤波器阶数选择。 cheb2ord -切比雪夫 (Chebyshev)II 型滤波器阶数 选择。 ellipord -椭圆 (Elliptic)滤波器阶数选择。 IIR 滤波器设计 butter -巴特沃斯 (Butterworth)滤波器设计 cheby1 -切比雪夫 (Chebyshev)I 型滤波器设计 cheby2 -切比雪夫 (Chebyshev)II 型滤波器设计 ellip -椭圆 (Elliptic)滤波器设计

3、maxflat -通用的巴特沃斯 (Butterworth)低通滤波器设计 yulewalk -Yule-Walker 滤波器设计 (直接数字滤波器设计法 ) 1. Butterworth 滤波器设计 Butterworth 滤波器是通带、阻带都单调衰减的滤波器。 (1)调用 buttord 函数确定巴特沃斯滤波器的阶数，格式为 N,Wc = buttord(Wp,Ws,Ap,As) 输入参数： Ap,As 为通带最大衰减和阻带最小衰减，以 dB 为单位。 Wp,Ws 为归一化通带截频和阻带截频， 0Wp,Ws1 。 输出参数： N 为滤波器的阶数； Wc 为截频， 0 Wc 1。 (2)调用

4、 butter 函数设计出巴特沃斯滤波器，格式为 iiNijjMjzazbzH101)(NNNNMMMM zazaza zbzbbb )1(111)1(11101 b,a = butter(N,Wc,options) 输入参数： N 和 Wc 是 buttord 函数返回的参数，含义见上。 Options=low, high, bandpass, stop, 分别对应 低通、高通、带通、带阻，默认情况下为低通或带通。 输出参数： b 和 a 为设计出的 IIR 数字滤波器 H(s)的分子多项式和分母多项式的系数矩阵。 2. Chebyshev II 型滤波器设计 Chebyshev II 型滤

5、波器为阻带纹波控制器： 在阻带呈现纹波特性。 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.信号 ,确定设计指标，实现各种 IIR 数字滤波器以实现以下信号处理。 （ 1）设计 IIR 低通滤波器，滤除 的成分。 （ 2）设计 IIR 高通滤波器，滤除 的成分。 （ 3）设计 IIR 带通滤波器，滤除 的成分。 （ 4）设

6、计 IIR 带阻滤波器，滤除 的成分。 要求利用 butterord 函数求解滤波器的阶数；利用 butter 函数设计各 IIR 数字滤波器；画出滤波器的幅度相应和相位响应；给出 IIR 数字滤波器的系统函数。 (1) clear; fsam=100; t0=1/fsam;t=6;k=0:t0:t; N=512; x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0); f=(-N/2:(N/2-1)/N*2*pi; N1,Wc=buttord(1/4,1/2,3,60); b,a=butter(N1,Wc,low); freqz(b,a); figure; axis(0

7、,1,-120,0); y=filter(b,a,x); Y=fftshift(fft(y,N); stem(f,abs(Y); )32cos()4cos(1 kkkx )32cos( k )4cos(1 k )32 k)4cos( k(2) clear; fsam=100; t0=1/fsam;t=6;k=0:t0:t; N=512; x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0); f=(-N/2:(N/2-1)/N*2*pi; N1,Wc=buttord(2/3,10/24,3,100); b,a=butter(N1,Wc,high); freqz(b,a);

8、 figure; axis(0,1,-120,0); y=filter(b,a,x); Y=fftshift(fft(y,N); stem(f,abs(Y); （ 3） clear; fsam=100; t0=1/fsam;t=6;k=0:t0:t; N=512; x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0); f=(-N/2:(N/2-1)/N*2*pi; N1,Wc=buttord(1/8 8/24,1/32 2/3,3,60); b,a=butter(N1,Wc,bandpass); freqz(b,a); figure; axis(0,1,-120,0);

9、 y=filter(b,a,x); Y=fftshift(fft(y,N); stem(f,abs(Y); (4) clear; fsam=100; t0=1/fsam;t=6;k=0:t0:t; N=512; x=1+cos(pi/4.*k/t0)+cos(2*pi/3.*k/t0); f=(-N/2:(N/2-1)/N*2*pi; N1,Wc=buttord(1/32 2/3,0.2 10/27,3,60); b,a=butter(N1,Wc,stop); freqz(b,a); figure; axis(0,1,-120,0); y=filter(b,a,x); Y=fftshift(f

10、ft(y,N); stem(f,abs(Y); 2.某带通滤波器的设计指标为 （ 1）试分别利用巴特沃斯、切比雪夫 型、切比雪夫 型和椭圆模拟滤波器，通过脉冲响应不变法设计该带通数字滤波器，画出其频率特性，比较设计结果。 （ 2）试分别利用巴特沃斯、切比雪夫 型、切比雪夫 型和椭圆模拟滤波器，通过双线性变换法设计该带通数字滤波器，画出其频率特性，比较设计结果。 （ 3）分析比较以上设计结果，有何结论？ (1).脉冲响应不变法 clear; omegas=0.2 0.72*pi; omegap=0.3 0.6*pi; Ap=1;As=42; Fs=1; ws=omegas*Fs; wp=omeg

11、ap*Fs; N,Wc=buttord(wp,ws,Ap,As,s); b,a=butter(N,Wc,bandpass,s); % N,Wc=cheb1ord(wp,ws,Ap,As,s); % b,a=cheby1(N,Ap,Wc,bandpass,s); % N,Wc=cheb2ord(wp,ws,Ap,As,s); % b,a=cheby2(N,As,Wc,bandpass,s); % dBAdBAradradradradspspps 42,1,)(72.0 ,)(6.0,)(3.0,)(2.02211 % N,Wc=ellipord(wp,ws,Ap,As,s); % b,a=ell

12、ip(N,Ap,As,Wc,bandpass,s); numd,dend=impinvar(b,a,Fs); w=linspace(0,pi,512); h=freqz(numd,dend,w); norm=max(abs(h); numd=numd/norm; plot(w/pi,20*log10(abs(h)/norm); grid on 右图从上到下分别是 buttord Cheby1 cheby2 ellipord 型滤波器的 频率特性 (脉冲响应不变法 ) (2)双线性变换法 clear; omegas=0.2 0.72*pi; omegap=0.3 0.6*pi; Ap=1;As=

13、42; Fs=1;T=1/Fs; wp=2*tan(omegap/2)/T; ws=2*tan(omegas/2)/T; % N,Wc=buttord(wp,ws,Ap,As,s); % b,a=butter(N,Wc,bandpass,s); % N,Wc=cheb1ord(wp,ws,Ap,As,s); % b,a=cheby1(N,Ap,Wc,bandpass,s); % % N,Wc=cheb2ord(wp,ws,Ap,As,s); % b,a=cheby2(N,As,Wc,bandpass,s); N,Wc=ellipord(wp,ws,Ap,As,s); b,a=ellip(N,A

14、p,As,Wc,bandpass,s); numd,dend=bilinear(b,a,Fs); w=linspace(0,pi,512); h=freqz(numd,dend,w); norm=max(abs(h); numd=numd/norm; plot(w/pi,20*log10(abs(h)/norm); grid on title(ellipord) （ 1）利用 DTMF 的产生原理合成数字“ 4”，合成公式为： ， 画出其时域波形和频谱图。 （ 2）根据图中的低通和高通滤波器，画出两个滤波器的幅频响应；并利用 对信号 进行滤波；利用 fft 函数求得输出的频响，绘出幅频响应曲线

15、。 （ 3）设计图中的第一组的四个带通滤波器，画出滤波器的幅频响应；并利用 对图中相应的输入信号进行滤波，在一张图上用不同颜色画出这四条输出曲线。比较其幅值， BP770Hz 滤波器的输出最大，因此确定该 DTMF 的行频为 770Hz。 (1) clear; k=0:204;N=205; x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k); stem(k,x); title(时域波形 ) figure; X=fftshift(fft(x); omega=(-(N-1)/2):(N-1)/2)/N*2*pi; stem(omega,abs(X); ti

16、tle(频谱图 ) )8000 12092sin()8000 7702sin(4 kkkx 204:0k ),( xabfiltery ),( xabfiltery 4x(2) clear; k=0:204;N=205; x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k); omega=(-(N-1)/2):(N-1)/2)/N*2*pi; %低通 wp=2*941/8000;ws=2*1200/8000;Ap=3;As=60; N,wc=buttord(wp,ws,Ap,As); b,a=butter(N,wc); freqz(b,a); figur

17、e %滤波输出 y=filter(b,a,x); Y=fftshift(fft(y); stem(omega,abs(Y) clear; k=0:204;N=205; x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k); omega=(-(N-1)/2):(N-1)/2)/N*2*pi; %高通 wp=2*1200/8000;ws=2*1000/8000;Ap=3;As=60; N,wc=buttord(wp,ws,Ap,As); b,a=butter(N,wc,high); freqz(b,a); figure %滤波输出 y=filter(b,a,

18、x); Y=fftshift(fft(y); stem(omega,abs(Y) (3) clear; k=0:204;N=205; x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k); omega=(-(N-1)/2):(N-1)/2)/N*2*pi; wp1=692 702./4000;ws1=687 707./4000; wp2=765 775./4000;ws2=760 780./4000; wp3=847 857./4000;ws3=842 862./4000; wp4=936 946./4000;ws4=931 951./4000; Ap=

19、3;As=30; N1,wc1=buttord(wp1,ws1,Ap,As); b1,a1=butter(N1,wc1,bandpass); N2,wc2=buttord(wp2,ws2,Ap,As); b2,a2=butter(N2,wc2,bandpass); N3,wc3=buttord(wp3,ws3,Ap,As); b3,a3=butter(N3,wc3,bandpass); N4,wc4=buttord(wp4,ws4,Ap,As); b4,a4=butter(N4,wc4,bandpass); figure(1) freqz(b1,a1); figure(2) freqz(b2,

20、a2); figure(3) freqz(b3,a3); figure(4) freqz(b4,a4); 以上分别是 4 个带通滤波器的幅频响应 . clear; k=0:204;N=205; x=sin(2*pi*770/8000.*k)+sin(2*pi*1209/8000.*k); omega=(-(N-1)/2):(N-1)/2)/N*2*pi; wp1=692 702./4000;ws1=687 707./4000; wp2=765 775./4000;ws2=760 780./4000; wp3=847 857./4000;ws3=842 862./4000; wp4=936 94

21、6./4000;ws4=931 951./4000; Ap=3;As=20; N1,wc1=buttord(wp1,ws1,Ap,As); b1,a1=butter(N1,wc1,bandpass); N2,wc2=buttord(wp2,ws2,Ap,As); b2,a2=butter(N2,wc2,bandpass); N3,wc3=buttord(wp3,ws3,Ap,As); b3,a3=butter(N3,wc3,bandpass); N4,wc4=buttord(wp4,ws4,Ap,As); b4,a4=butter(N4,wc4,bandpass); y1=filter(b1,a1,x); Y1=abs(fftshift(fft(y1); y2=filter(b2,a2,x); Y2=abs(fftshift(fft(y2); y3=filter(b3,a3,x); Y3=abs(fftshift(fft(y3); y4=filter(b4,a4,x); Y4=abs(fftshift(fft(y4); plot(omega,Y1,-g,omega,Y2,-r,omega,Y3,-b,omega,Y4,-k);