数字信号处理实验七

1、实 验 报 告 实验名称:IIR数字滤波器设计及应用 课程名称: 数字信号处理 院 系 部:电气与电子工程学院 专业班级：信息1002学生姓名：王萌 学 号: 1101200219同 组 人: 实验台号: 指导教师：范杰清 成 绩： 实验日期: 华北电力大学（北京）一、实验目的及要求：加深理解IIR数字滤波器的特性，掌握IIR数字滤波器的设计原理与设计方法，以及IIR数字滤波器的应用。二、实验原理N阶IIR数字滤波器的系统函数为: IIR数字滤波器的设计主要通过成熟的模拟滤波器设计方法来实现：将数字滤波器设计指标转换为模拟滤波器设计指标，设计出相应的模拟滤波器H(s)，再经过脉冲响应不变法或双

2、线性变换法得到所需的IIR数字滤波器H(z)。IIR数字滤波器设计的重要环节是模拟原型低通滤波器的设计，主要包括Butterworth、Chebyshev和椭圆等滤波器。MATLAB 信号处理工具箱中提供了IIR滤波器设计的函数。n IIR 滤波器阶数选择buttord -巴特沃斯(Butterworth)滤波器阶数选择。cheb1ord -切比雪夫(Chebyshev)I 型滤波器阶数选择。cheb2ord -切比雪夫(Chebyshev)II 型滤波器阶数选择。ellipord -椭圆(Elliptic)滤波器阶数选择。n IIR 滤波器设计butter -巴特沃斯(Butterworth

3、)滤波器设计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为通带最大衰减和阻带最小衰减，以

4、dB为单位。Wp,Ws为归一化通带截频和阻带截频，0<Wp,Ws<1 。输出参数：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型

5、滤波器为通带纹波控制器：在通带呈现纹波特性，在阻带单调衰减。N,Wc = cheb1ord(Wp, Ws, Ap, As)b,a = cheby1(N,Ap,Wc,options)参数含义与butter中参数一致。3. Chebyshev II 型滤波器设计Chebyshev II 型滤波器为阻带纹波控制器：在阻带呈现纹波特性。N,Wc = cheb2ord(Wp, Ws, Ap, As)b,a = cheby2(N,As,Wc,options)4. 椭圆滤波器设计椭圆滤波器在通阻带都呈现纹波特性。N,Wc = ellipord(Wp,Ws,Ap,As)b,a = ellip(N,Ap,As,

6、Wc,options)三、实验例题1.设计一个数字低通滤波器，要求在00.2p(rad)内衰耗不大于3dB，在0.6(rad)内衰耗不小于60dB。若模拟低通原型采用Butterworth滤波器：N,Wc=buttord(0.2,0.6,3,60);b,a=butter(N,Wc);freqz(b,a);axis(0,1,-120,0);grid on 2.设计一个数字带通滤波器，在100Hz 200Hz通带内纹波不大于3dB，通带两边各50 Hz外是阻带，衰耗不小于40dB。抽样频率为1000 Hz 。Wp =100 200/500;Ws =100-50 200+50/500; Rp = 3

7、; Rs = 40;N,Wn=cheb1ord(Wp,Ws,Rp,Rs)b,a =cheby1(N,Rp,Wn);freqz(b,a,512,1000); title('Chebyshev Type I Bandpass Filter')axis(0,500,-80,0); grid on 五、实验内容1.信号 ,确定设计指标，实现各种IIR数字滤波器以实现以下信号处理。要求利用butterord函数求解滤波器的阶数；利用butter函数设计各IIR数字滤波器；画出滤波器的幅度相应和相位响应；给出IIR数字滤波器的系统函数。（1）设计IIR低通滤波器，滤除的成分。N,Wc=bu

8、ttord(0.25,2/3,3,60);b,a=butter(N,Wc);freqz(b,a,512,1000);title('n=5 Butterworth Lowpass Filter');N=0:1:100;x=1+cos(pi*N/4)+cos(2*N*pi/3);figure(2),subplot(1,2,1);stem(N,abs(fft(x);title('原信号频谱');subplot(1,2,2);stem(N,angle(fft(x);y=filter(b,a,x);figure(3);stem(N,y);title('滤波后信号时

9、域波形');figure(4),subplot(1,2,1);stem(N,abs(fft(y);title('滤波后信号频谱');axis(0,100,0,60); subplot(1,2,2);stem(N,angle(fft(y); （2）设计IIR高通滤波器，滤除的成分。N,Wc=buttord(0.66,0.25,3,60);b,a=butter(N,Wc,'high');freqz(b,a,512,1000);title('n=5 Butterworth highpass Filter');N=0:1:100;x=1+cos(

10、pi*N/4)+cos(2*N*pi/3);figure(2),subplot(1,2,1);stem(N,abs(fft(x);title('原信号频谱');subplot(1,2,2);stem(N,angle(fft(x);y=filter(b,a,x);figure(3);stem(N,y);title('滤波后信号时域波形');figure(4),subplot(1,2,1);stem(N,abs(fft(y);title('滤波后信号频谱');axis(0,100,0,60); subplot(1,2,2);stem(N,angle(

11、fft(y);（3） 设计IIR带通滤波器，滤除的成分。Wp =0.25 0.6; Ws =0.25-0.2 0.6+0.2; Rp = 3; Rs = 60;N,Wc=buttord(Wp,Ws,Rp,Rs)b,a=butter(N,Wc,'bandpass')freqz(b,a,512,1000);title('n=5 Butterworth bandpass Filter');N=0:1:100;x=1+cos(pi*N/4)+cos(2*N*pi/3);figure(2),subplot(1,2,1);stem(N,abs(fft(x);title(&#

12、39;原信号频谱');subplot(1,2,2);stem(N,angle(fft(x);y=filter(b,a,x);figure(3);stem(N,y);title('滤波后信号时域波形');figure(4),subplot(1,2,1);stem(N,abs(fft(y);title('滤波后信号频谱');axis(0,100,0,60); subplot(1,2,2);stem(N,angle(fft(y);（4） 设计IIR带阻滤波器，滤除的成分。Ws =0.25 0.5; Wp =0.25-0.1 0.5+0.1; Rp = 3; R

13、s = 60;N,Wc=buttord(Wp,Ws,Rp,Rs)b,a=butter(N,Wc,'stop')freqz(b,a,512,1000);title('n=5 Butterworth stop Filter');N=0:1:100;x=1+cos(pi*N/4)+cos(2*N*pi/3);figure(2),subplot(1,2,1);stem(N,abs(fft(x);title('原信号频谱');subplot(1,2,2);stem(N,angle(fft(x);y=filter(b,a,x);figure(3);stem(

14、N,y);title('滤波后信号时域波形');figure(4),subplot(1,2,1);stem(N,abs(fft(y);title('滤波后信号频谱');axis(0,100,0,60); subplot(1,2,2);stem(N,angle(fft(y);2.某带通滤波器的设计指标为 （1）试分别利用巴特沃斯、切比雪夫型、切比雪夫型和椭圆模拟滤波器，通过脉冲响应不变法设计该带通数字滤波器，画出其频率特性，比较设计结果。脉冲响应不变法设计 巴特沃斯 带通数字滤波器Wp =0.3 0.6; Ws =0.2 0.72; Ap = 1; As = 42

15、; Fs=1;N,Wc=buttord(Wp,Ws,Ap,As);b,a=butter(N,Wc,'bandpass');h=freqz(b,a,w);axis(0,1,-120,0);numd,dend=impinvar(b,a,Fs);w=linspace(0,pi,512);norm=max(abs(h);a=a/norm;plot(w/pi,20*log10(abs(h)/norm); title('脉冲响应不变法 巴特沃斯 带通数字滤波器') 脉冲响应不变法设计 切比雪夫1型 带通数字滤波器Wp =0.3 0.6; Ws =0.2 0.72; Ap =

16、 1; As = 42; Fs=1;N,Wc=cheb1ord(Wp,Ws,Ap,As);b,a=cheby1(N,As,Wc,'bandpass');w=linspace(0,pi,512); h=freqz(b,a,w);axis(0,1,-120,0);numd,dend=impinvar(b,a,Fs);w=linspace(0,pi,512);norm=max(abs(h); a=a/norm; plot(w/pi,20*log10(abs(h)/norm);title('脉冲响应不变法 切比雪夫1型 带通数字滤波器');脉冲响应不变法设计 切比雪夫2

17、型 带通数字滤波器Wp =0.3 0.6; Ws =0.2 0.72; Ap = 1; As = 42; Fs=1;N,Wc=cheb2ord(Wp,Ws,Ap,As);b,a=cheby2(N,As,Wc,'bandpass');w=linspace(0,pi,512);h=freqz(b,a,w); axis(0,1,-120,0);numd,dend=impinvar(b,a,Fs);norm=max(abs(h); a=a/norm;plot(w/pi,20*log10(abs(h)/norm);title('脉冲响应不变法 切比雪夫2型 带通数字滤波器'

18、;);脉冲响应不变法设计 椭圆 型带通数字滤波器Wp =0.3 0.6; Ws =0.2 0.72; Ap = 1; As = 42; Fs=1;N,Wc=ellipord(Wp,Ws,Ap,As);b,a=ellip(N,Ap,As,Wc,'bandpass');h=freqz(b,a,w); axis(0,1,-120,0);numd,dend=impinvar(b,a,Fs);w=linspace(0,pi,512);norm=max(abs(h); a=a/norm;plot(w/pi,20*log10(abs(h)/norm);title('脉冲响应不变法 椭

19、圆 带通数字滤波器');（2） 试分别利用巴特沃斯、切比雪夫型、切比雪夫型和椭圆模拟滤波器，通过双线性变换法设计该带通数字滤波器，画出其频率特性，比较设计结果。双线性变换法设计巴特沃斯 带通数字滤波器WWp =0.3 0.6; WWs =0.2 0.72; Wp=tan(WWp/2);Ws=tan(WWs/2);Ap = 1; As = 42; Fs=1;N,Wc=buttord(Wp,Ws,Ap,As);b,a=butter(N,Wc,'bandpass');h=freqz(b,a,w); axis(0,1,-120,0);numd,dend=impinvar(b,a

20、,Fs);w=linspace(0,pi,512);norm=max(abs(h); a=a/norm;plot(w/pi,20*log10(abs(h)/norm);title('双线性变换法 巴特沃斯 带通数字滤波器')双线性变换法设计 切比雪夫1型 带通数字滤波器WWp =0.3 0.6; WWs =0.2 0.72; Wp=tan(WWp/2);Ws=tan(WWs/2);Ap = 1; As = 42; Fs=1;N,Wc=cheb1ord(Wp,Ws,Ap,As);b,a=cheby1(N,As,Wc,'bandpass');w=linspace(0

21、,pi,512);h=freqz(b,a,w);axis(0,1,-120,0);numd,dend=impinvar(b,a,Fs);w=linspace(0,pi,512);norm=max(abs(h); a=a/norm;plot(w/pi,20*log10(abs(h)/norm);title('双线性变换法 切比雪夫1型 带通数字滤波器');双线性变换法设计 切比雪夫2型 带通数字滤波器WWp =0.3 0.6; WWs =0.2 0.72; Wp=tan(WWp/2);Ws=tan(WWs/2);Ap = 1; As = 42; Fs=1;N,Wc=cheb2or

22、d(Wp,Ws,Ap,As);b,a=cheby2(N,As,Wc,'bandpass');w=linspace(0,pi,512);h=freqz(b,a,w);axis(0,1,-120,0);numd,dend=impinvar(b,a,Fs);w=linspace(0,pi,512);norm=max(abs(h); a=a/norm;plot(w/pi,20*log10(abs(h)/norm);title('双线性变换法 切比雪夫2型 带通数字滤波器');双线性变换法设计 椭圆 型带通数字滤波器WWp =0.3 0.6; WWs =0.2 0.72;

23、 Wp=tan(WWp/2);Ws=tan(WWs/2);Ap = 1; As = 42; Fs=1;N,Wc=ellipord(Wp,Ws,Ap,As);b,a=ellip(N,Ap,As,Wc,'bandpass');h=freqz(b,a,w); axis(0,1,-120,0);numd,dend=impinvar(b,a,Fs);w=linspace(0,pi,512);norm=max(abs(h); a=a/norm;plot(w/pi,20*log10(abs(h)/norm);title('双线性变换法 椭圆 带通数字滤波器');（3）分析比较

24、以上设计结果，有何结论？双线性变换法设计的滤波器数字频率和模拟频率之间是非线性关系3. 双音多频( DTMF)信号的产生及检测在现代通信系统中有着广泛的应用，家用电话、移动电话以及公共程控交换机(PBX)都采用DTMF信号发送和接收电话拨号号码。DTMF信号由低频组和高频组两组频率信号构成，高频为列，低频为行。按键电话上每个按键都由对应的两个频率组成。当按下某个键时，所得到的按键信号是由相应两个频率的正弦信号相加而成： xk=sin(k)+sin(k)， ，。DTMF信号可以用IIR数字滤波器来进行解码。如图所示接收到的音频信号通过低通和高通滤波器。由于DTMF信号行频低于列频，所以低通和高通

25、可以分离出行和列信息。再通过带通滤波器，对其输出的振幅大小进行判断，可以确定具体的行频和列频。如图中所示的滤波器响应，将得到按键“4”。（1）利用DTMF的产生原理合成数字“4”，合成公式为： ，画出其时域波形和频谱图。T0=204; N=205; T=1; k=0:T: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

26、,abs(fftshift(Xm);xlabel('f(Hz)');ylabel('magnitude');title('幅度谱');（2） 根据图中的低通和高通滤波器，画出两个滤波器的幅频响应；并利用对信号进行滤波；利用fft函数求得输出的频响，绘出幅频响应曲线。N,Wc=buttord(0.2,0.3,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:T:T0;x=sin

27、(2/8000)*770*pi*k)+sin(2/8000)*1209*pi*k); y=filter(b,a,x);figure(2)subplot(2,1,1);stem(k,y);title('对信号进行低通滤滤波')ym=fft(y,N)/N; subplot(2,1,2);stem(k,abs(fftshift(ym);title('低通滤滤波后频谱图')N,Wc=buttord(0.3,0.2,3,60); b,a=butter(N,Wc,'high');freqz(b,a); axis(0,1,-120,0);grid on;tit

28、le('巴特沃斯高通数字滤波器')T0=204; N=205; T=1; k=0:T:T0;x=sin(2/8000)*770*pi*k)+sin(2/8000)*1209*pi*k); y=filter(b,a,x);figure(2);subplot(2,1,1);stem(k,y);title('对信号进行高通滤滤波')ym=fft(y,N)/N; subplot(2,1,2);stem(k,abs(fftshift(ym);title('高通滤滤波后频谱图')（3） 设计图中的第一组的四个带通滤波器，画出滤波器的幅频响应；并利用对图中相应

29、的输入信号进行滤波，在一张图上用不同颜色画出这四条输出曲线。比较其幅值，BP770Hz滤波器的输出最大，因此确定该DTMF的行频为770Hz。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('Chebyshev Type I Bandpass1 Filter'); axis(0,4000,-80,0); grid onWp2 =7

30、50 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('Chebyshev Type I Bandpass2 Filter'); axis(0,4000,-80,0); grid onWp3 =830 870/4000; Ws3 =830-20 870+20/4000; Rp3 = 3; Rs3 = 40; N3,Wn3=cheb

31、1ord(Wp3,Ws3,Rp3,Rs3)b3,a3 =cheby1(N3,Rp3,Wn3); figure; freqz(b3,a3,512,8000); title('Chebyshev Type I Bandpass3 Filter'); axis(0,4000,-80,0); 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

32、,512,8000); title('Chebyshev Type I Bandpass4 Filter'); axis(0,4000,-80,0); grid on; k=0:1:204;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,'r-',k,y3,'y-',k,y4,'g-'

33、;); title('输出曲线') ; legend('697HZ', '770HZ','852HZ','941HZ');（4） 设计图中的第二组的三个带通滤波器，画出滤波器的幅频响应；并利用 对图中相应的输入信号进行滤波，在一张图上用不同颜色画出这三条输出曲线。比较其幅值，BP1209Hz滤波器的输出最大，因此确定该DTMF的列频为1209Hz。从（4）和（5）可以解出该DTMF信号为按键“4”。Wp1 =1180 1220/4000; Ws1 =1180-30 1220+30/4000; Rp1 = 3; Rs

34、1 = 40; N1,Wn1=cheb1ord(Wp1,Ws1,Rp1,Rs1)b1,a1 =cheby1(N1,Rp1,Wn1); freqz(b1,a1,512,8000); title('Chebyshev Type I Bandpass1 Filter'); axis(0,4000,-80,0); 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); fig

35、ure; freqz(b2,a2,512,8000); title('Chebyshev Type I Bandpass2 Filter'); axis(0,4000,-80,0); 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('Chebyshev Type I Bandpass

36、3 Filter'); axis(0,4000,-80,0); grid on;k=0:1:204; 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,'r-',k,y3,'g-'); title('输出曲线') ; legend('1209HZ','1336HZ','1477HZ')

37、;（5） 生成09其他数字的DTMF信号，按照上述方法进行解码，验证该方法的正确性。k=0:1:204; 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

38、=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 =

39、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(N

40、3,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 13

41、50+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

42、,x0); figure ; subplot(2,1,1); plot(k,Y01,k,Y02,'r-',k,Y03,'y-',k,Y04,'g-'); title('输出曲线0');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,'

43、;r-',k,y03,'g-'); 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,'r-',k,Y13,'y-',k,Y14,'g-'); title('输出曲线1') ; legend(&

44、#39;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,'r-',k,y13,'g-'); legend('1209HZ', '1336HZ','1477HZ');Y21=filter(B1,A1,x2);Y22=filter(B2,A2,x2)

45、;Y23=filter(B3,A3,x2);Y24=filter(B4,A4,x2);figure; subplot(2,1,1); plot(k,Y21,k,Y22,'r-',k,Y23,'y-',k,Y24,'g-'); title('输出曲线2') ;legend('697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y21=filter(b1,a1,x2); y22=filter(b2,a2,x2); y23=f

46、ilter(b3,a3,x2); plot(k,y21,k,y22,'r-',k,y23,'g-'); legend('1209HZ', '1336HZ','1477HZ');Y31=filter(B1,A1,x3);Y32=filter(B2,A2,x3);Y33=filter(B3,A3,x3);Y34=filter(B4,A4,x3);figure;subplot(2,1,1); plot(k,Y31,k,Y32,'r-',k,Y33,'y-',k,Y34,'g-

47、9;); figure; subplot(2,1,1); plot(k,Y21,k,Y22,'r-',k,Y23,'y-',k,Y24,'g-'); title('输出曲线3') ; legend('697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y31=filter(b1,a1,x3); y32=filter(b2,a2,x3); y33=filter(b3,a3,x3); plot(k,y31,k,y32,

48、9;r-',k,y33,'g-'); legend('1209HZ', '1336HZ','1477HZ');Y41=filter(B1,A1,x4);Y42=filter(B2,A2,x4);Y43=filter(B3,A3,x4);Y44=filter(B4,A4,x4);figure;subplot(2,1,1); plot(k,Y41,k,Y42,'r-',k,Y43,'y-',k,Y44,'g-'); title('输出曲线4') ;legend(&

49、#39;697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y41=filter(b1,a1,x4); y42=filter(b2,a2,x4); y43=filter(b3,a3,x4); plot(k,y41,k,y42,'r-',k,y43,'g-'); legend('1209HZ', '1336HZ','1477HZ');Y51=filter(B1,A1,x5);Y52=filter(B2,A2,x5)

50、;Y53=filter(B3,A3,x5);Y54=filter(B4,A4,x5);figure;subplot(2,1,1); plot(k,Y51,k,Y52,'r-',k,Y53,'y-',k,Y54,'g-'); title('输出曲线5') ; legend('697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y51=filter(b1,a1,x5); y52=filter(b2,a2,x5); y53=f

51、ilter(b3,a3,x5); plot(k,y51,k,y52,'r-',k,y53,'g-'); legend('1209HZ', '1336HZ','1477HZ');Y61=filter(B1,A1,x6);Y62=filter(B2,A2,x6);Y63=filter(B3,A3,x6);Y64=filter(B4,A4,x6);figure;subplot(2,1,1); plot(k,Y61,k,Y62,'r-',k,Y63,'y-',k,Y64,'g-

52、9;); title('输出曲线6') ; legend('697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y61=filter(b1,a1,x6); y62=filter(b2,a2,x6); y63=filter(b3,a3,x6); plot(k,y61,k,y62,'r-',k,y63,'g-'); legend('1209HZ', '1336HZ','1477HZ');Y71

53、=filter(B1,A1,x7);Y72=filter(B2,A2,x7);Y73=filter(B3,A3,x7);Y74=filter(B4,A4,x7);figure;subplot(2,1,1); plot(k,Y71,k,Y72,'r-',k,Y73,'y-',k,Y74,'g-'); title('输出曲线7') ; legend('697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y71=filter(

54、b1,a1,x7); y72=filter(b2,a2,x7); y73=filter(b3,a3,x7); plot(k,y71,k,y72,'r-',k,y73,'g-'); legend('1209HZ', '1336HZ','1477HZ');Y81=filter(B1,A1,x8);Y82=filter(B2,A2,x8);Y83=filter(B3,A3,x8);Y84=filter(B4,A4,x8);figure;subplot(2,1,1); plot(k,Y81,k,Y82,'r-

55、9;,k,Y83,'y-',k,Y84,'g-'); title('输出曲线8'); legend('697HZ', '770HZ','852HZ','941HZ'); subplot(2,1,2); y81=filter(b1,a1,x8); y82=filter(b2,a2,x8); y83=filter(b3,a3,x8); plot(k,y81,k,y82,'r-',k,y83,'g-'); legend('1209HZ', &

56、#39;1336HZ','1477HZ');Y91=filter(B1,A1,x9);Y92=filter(B2,A2,x9);Y93=filter(B3,A3,x9);Y94=filter(B4,A4,x9);figure;subplot(2,1,1); plot(k,Y91,k,Y92,'r-',k,Y93,'y-',k,Y94,'g-'); title('输出曲线9') ; legend('697HZ', '770HZ','852HZ','941

57、HZ'); subplot(2,1,2); y91=filter(b1,a1,x9); y92=filter(b2,a2,x9); y93=filter(b3,a3,x9); plot(k,y91,k,y92,'r-',k,y93,'g-'); legend('1209HZ', '1336HZ','1477HZ');（6） 对于图中的第一组带通滤波器的4个输出，分别用fft函数求其频响，画出幅频响应曲线。同理在一张图上用不同颜色画出第二组带通滤波器的3个输出的幅频响应曲线。能否用这种方法解出DTMF信号？T

58、0=204; N=205; 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)