《数字信号处理》课程研究性学习报告IIR和FIR滤波器设计专题研讨

1、数字信号处理课程研究性学习报告姓名 学号 同组成员 指导教师 时间 2011/5/10 iir和fir滤波器设计专题研讨【目的】(1) 掌握根据滤波器指标设计iir和fir数字滤波器的原理和方法。(2) 熟悉通过iir和fir数字滤波器进行实际系统设计的方法。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】 1 设计一个数字滤波器，在频率低于的范围内，低通幅度特性为常数，并且不低于0.75db。在频率和之间，阻带衰减至少为20db。(1)试求满足这些条件的最低阶butterworth滤波器。(2)试求满足这些条件的最低阶chebyshev i滤波器。(3)自主

2、选择一段带限信号，通过所设计的（1）、（2）两种滤波器，比较各自的输入和输出信号。讨论两种滤波器在结构和性能上的差异。【题目分析】本题讨论模拟滤波器和数字滤波器的设计。结合课本分析各类滤波器的特性。【iir模拟滤波器设计的基本方法】在设计iir模拟滤波器的时候，1） 先将待设计的模拟滤波器技术指标转换为原型模拟低通滤波器的指标；2） 然后设计原型模拟低通滤波器（有butterworth、chebyshev和椭圆等设计方法）；3） 再通过频率变换将模拟低通滤波器转换为所需的滤波器。【仿真结果】【结果分析】butterworth滤波器:从图上可以看出，所设计的滤波器满足指标，通带内为恒定常数0db

3、，在过渡带内单调递减，在pi附近急剧下降，阻带衰减有裕量。chebyshev i滤波器：从图上可以看出，所设计的滤波器满足指标，通带波动，过渡带内单调递减，相对于butterworth滤波器的阻带衰减比较均匀。【自主学习内容】1） 用butterworth、chebyshev和椭圆等设计方法设计模拟低通滤波器；2） 如何用matlab实现上述模拟滤波器的实现；3） 通过频率变换将模拟低通滤波器转换为所需的滤波器；4） 用单位脉冲响应法和双线性变换法实现模拟滤波器向数字滤波器的转换。【阅读文献】数字信号处理（第二版） 陈后金主编 高等教育出版社【发现问题】 (专题研讨或相关知识点学习中发现的问题

4、)：无【问题探究】单位脉冲响应法与双线性变换法的比较：（使用butterworth的设计方法讨论）仿真程序wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=buttord(wp,ws,ap,as,s);num_a,den_a=butter(n,wc,s);num_d1,den_d1=impinvar(num_a,den_a,fs);num_d2,den_d2=bilinear(num_a,den_a,fs);w=linspace(0,pi,1024);h1=freqz(num_d1,den_d1,w

5、);norm1=max(abs(h1);num_d1=num_d1/norm1;h2=freqz(num_d2,den_d2,w);norm2=max(abs(h2);num_d2=num_d2/norm2;plot(w/pi,20*log10(abs(h1/norm1),blue);hold on;plot(w/pi,20*log10(abs(h2/norm2),green);xlabel(normalized frequency);ylabel(gain,db);仿真结果结果分析在本次设计中，通带内，用单位脉冲响应法和双线性变换区别并不明显，相比于单位脉冲响应法，双线性变换法，通带截频处地

6、衰减略大，在过渡带和阻带则出现明显的差别，即双线性变换法的衰减明显大于单位脉冲响应法。【仿真程序】butterworth滤波器:wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=buttord(wp,ws,ap,as,s);num_a,den_a=butter(n,wc,s);num_d,den_d=impinvar(num_a,den_a,fs);w=linspace(0,pi,1024);h=freqz(num_d,den_d,w);norm=max(abs(h);num_d=num_d/nor

7、m;plot(w/pi,20*log10(abs(h/norm);xlabel(normalized frequency);ylabel(gain,db);chebyshev i滤波器：wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=cheb1ord(wp,ws,ap,as,s);num_a,den_a=cheby1(n,ap,wc,s);num_d,den_d=impinvar(num_a,den_a,fs);w=linspace(0,pi,1024);h=freqz(num_d,den_d,w

8、);norm=max(abs(h);num_d=num_d/norm;plot(w/pi,20*log10(abs(h/norm);xlabel(normalized frequency);ylabel(gain,db);滤波：1）butterworth滤波器wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=buttord(wp,ws,ap,as,s);num_a,den_a=butter(n,wc,s);num_d,den_d=impinvar(num_a,den_a,fs); w=linspac

9、e(0,pi,1024);h=freqz(num_d,den_d,w);x=w;s=h.*w;plot(w,abs(s);2）chebyshev i滤波器wp=0.2613*pi; ws=0.4018*pi; ap=0.75; as=20;fs=1; wp=wp*fs; ws=ws*fs;n,wc=cheb1ord(wp,ws,ap,as,s);num_a,den_a=cheby1(n,ap,wc,s);num_d,den_d=impinvar(num_a,den_a,fs);w=linspace(0,pi,1024);h=freqz(num_d,den_d,w);x=w;s=h.*w;plo

10、t(w,abs(s);2 分别用hamming窗，blackman窗和kaiser窗设计，满足下列指标的fir低通滤波器：(1)画出所设计滤波器的幅频响应。(2)比较这三种窗的设计结果。【题目分析】本题讨论窗函数法设计数字fir滤波器。结合课本分析不同窗函数法的设计结果。【fir模拟滤波器设计的基本方法】1、 窗函数法1） 根据待逼近的理想滤波器的频率响应，确定滤波器的类型；2） 确定理想滤波器的幅度函数和相位函数，由idtft求出理想滤波器的单位脉冲响应；3） 将上述所求的无限长单位脉冲响应加窗截断得到有限长的序列。2、 频率取样法1） 根据待逼近的理想滤波器的频率响应，确定滤波器的类型；2

11、） 确定理想滤波器的幅度函数、相位函数和频率响应上的取样点，由idft得到所需的滤波器有限长的单位脉冲响应。【仿真结果】【结果分析】从上面的仿真可以看出，1） 分别用hamming窗、blackman窗、kaiser窗设计的fir低通滤波器在通带内没有波动；2） 通带衰减和阻带衰减，kaiser窗最小，hamming窗次之，blackman窗最大；3） 三种窗设计的阶数由kaiser窗、hamming窗、blackman窗依次增加。【自主学习内容】1) 窗函数设计线性相位fir滤波器的思想和步骤；2) 各种窗函数的具体内容、区别和使用；3) 频率频率取样法设计线性相位fir滤波器的思想及步骤；

12、4) 用matlab实现上述两种方法的设计方法。【阅读文献】数字信号处理（第二版） 陈后金主编 高等教育出版社【发现问题】 (专题研讨或相关知识点学习中发现的问题)：kaiser窗是一种可调窗，如何调节阻带衰减？【问题探究】kaiser窗rs参数的调节对滤波器设计的影响（程序与）仿真结果 结果分析1） rs影响阻带衰减，rs越小，阻带衰减越大；2） rs影响滤波器的阶数，rs越大，阶数越高；3） rs影响滤波器的阻带波动范围。【仿真程序】hamming窗:（n=35）wp=0.4*pi;ws=0.6*pi;ap=0.5;as=45;n=ceil(7*pi/(ws-wp);n=mod(n+1,2

13、)+n;m=n-1;fprintf(n=%.0fn,n);w=hamming(n);wc=(wp+ws)/2;k=0:m;hd=(wc/pi)*sinc(wc*(k-0.5*m)/pi);h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb);blackman窗:（n=59）wp=0.4*pi;ws=0.6*pi;ap=0.5;as=45;n=ceil(11.4*pi/(ws-wp);n=mod(n+1,2)+n;m=n-1;fprintf(n=%.0f

14、n,n);w=blackman(n);wc=(wp+ws)/2;k=0:m;hd=(wc/pi)*sinc(wc*(k-0.5*m)/pi);h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb);kaiser窗:（n=24）rs=0.01;f=0.4 0.6;a=1,0;dev=rs*ones(1,length(a);m,wc,beta,ftype=kaiserord(f,a,dev);h=fir1(m,wc,ftype,kaiser(m+1,beta

15、);omega=linspace(0,pi,512);mag=freqz(h,1,omega);plot(omega/pi,20*log10(abs(mag);三种窗：rs=0.01;f=0.4 0.6;a=1,0;dev=rs*ones(1,length(a);m1,wc1,beta,ftype=kaiserord(f,a,dev);h1=fir1(m1,wc1,ftype,kaiser(m1+1,beta);omega=linspace(0,pi,512);mag1=freqz(h1,1,omega);plot(omega/pi,20*log10(abs(mag1),blue);hold

16、on;wp=0.4*pi;ws=0.6*pi;ap=0.5;as=45;n2=ceil(7*pi/(ws-wp);n2=mod(n2+1,2)+n;m2=n2-1;fprintf(n2=%.0fn,n2);w2=hamming(n2);wc2=(wp+ws)/2;k=0:m2;hd2=(wc2/pi)*sinc(wc2*(k-0.5*m2)/pi);h2=hd2.*w2;mag2=freqz(h2,1,omega);magdb2=20*log10(abs(mag2);plot(omega/pi,magdb2,green);hold on;wp=0.4*pi;ws=0.6*pi;ap=0.5;a

17、s=45;n3=ceil(11.4*pi/(ws-wp);n3=mod(n3+1,2)+n3;m3=n3-1;fprintf(n3=%.0fn,n3);w3=blackman(n3);wc3=(wp+ws)/2;k=0:m3;hd3=(wc3/pi)*sinc(wc3*(k-0.5*m3)/pi);h3=hd3.*w3;mag3=freqz(h3,1,omega);magdb3=20*log10(abs(mag3);plot(omega/pi,magdb3,red);3附件给出了一段含有噪声的音频信号。（1）分析该信号的频谱特点。（2）通过脉冲响应不变法设计一个滤波器对其进行处理，得到有用信息

18、，自主确定各项指标。（3）使用双线性变换法重做（2）。（4）试定量比较上述两种滤波器的各项性能，画出能说明性能差异的相关图形，对比并解释。（5）通过不同窗函数法设计fir滤波器对其进行处理，比较不同窗函数的滤波效果，给出分析和结论。（选做）(6) 请尝试采用其它的音频信号，混入不同的噪声，利用所学的滤波方法进行分析，会得到什么样的效果？【题目分析】本题讨论用iir和fir数字滤波器进行实际系统设计的方法。【仿真结果】(1)【结果分析】对几种滤波器得到的音频信号进行对比。【自主学习内容】这条指令主要功能是读取音频文件，该音频文件的格式规定为microsoft wave格式，也就是以.wav为后缀

19、的音频文件。有以下几种使用方法。1、 y = wavread(file);该语句读取文件名为file的音频文件，并将音频数据存储在y中，返回给用户使用。音频数据的值域为-1,12、y,fs,nbits = wavread(file)同上，file指定音频文件名，y是音频数据。另外该语句还能够返回该音频文件的采样率fs,和采样位数nbits3、. = wavread(file,n)该语句用来读取文件名为file的音频文件中的前n个数据。如果有多个通道，每个通道的数据都会被读取出来。.表示随用户自己定义。如果用户在这里只写一个y，该语句就只会返回音频数据，如果用户写成y,fs,nbits,该语句就

20、不但会返回音频数据y,还会返回采样率fs，和采样位数nbits4、.=wavread(file,n1,n2)该语句用来读取文件名为file的音频文件中的第n1个到第n2个之间的数据。5、siz = wavread(file,size)该语句返回文件名为file的音频文件的数据大小 siz = samples channels6、y,fs,nbits,opts=wavread(.) 该语句返回的参数opts包含了音频文件的额外信息。【阅读文献】数字信号处理（第二版） 陈后金主编 高等教育出版社【发现问题】 (专题研讨或相关知识点学习中发现的问题)：当有裕量出现时，选择不同的窗对滤波器选择有阶数的

21、影响。【问题探究】在iir滤波器设计过程中，由于利用指标数据的不同，造成裕量的出现。讨论利用不同指标出现裕量对滤波器性能的影响，以及如何有效地利用它？在fir滤波器设计过程中，由于所选用窗函数的不同，导致对信号滤波的效果不同。本题的语音信号加入的是双频噪声，讨论如果对信号加入其它形式的噪声，采用何种滤波器更合适。【仿真程序】（1）含有噪声音频分析:x,fs,bits=wavread(hello.wav); x1=x(1:1:end);x1=fft(x1); plot(abs(x1);title(含有噪声的音频信号);仿真结果结果分析读取音频文件后，得知抽样频率为fs=44100hz,共取n=5

22、0688个点。经过将图放大，在8800hz8950hz，1005hz1025hz两处范围内有噪声信号，发现所需要设计的滤波器参数应满足wp6000,ws8800,才能更可能将噪声去除。(2) 利用脉冲响应不变法所设计的低通滤波器: wp=6336;ws=8448;ap=1;as=25;（让角频率可以取整数）x,fs,bits=wavread(hello.wav); x1=x(1:1:end);x1=fft(x1); w=linspace(0,2*pi,50688);subplot(3,1,1);plot(w/pi,abs(x1);title(含有噪声的音频信号频谱图);np=6336;ns=8

23、448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;subplot(3,1,2);plot(w/pi,20*log10(abs(h)/norm); title(利用脉冲响应不变法

24、所设计的滤波器幅度分布);x=x1.*(h);subplot(3,1,3);plot(w/pi,abs(x);title(滤除噪声后的音频信号频谱图);仿真结果结果分析要使噪声消除，则滤波器在噪声频率处需要完成衰减，此音频信号中，需要完成的衰减约为49.54db，(3) 利用双线性变换法所设计的低通滤波器: wp=6336;ws=8448;ap=1;as=25;（让角频率可以取整数）x,fs,bits=wavread(hello.wav); x1=x(1:1:end);x1=fft(x1); w=linspace(0,2*pi,50688);subplot(3,1,1);plot(w/pi,a

25、bs(x1);title(含有噪声的音频信号频谱图);np=6336;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;subplot(3,1,2);plot(w/pi

26、,20*log10(abs(h)/norm); title(利用双线性变化法所设计的滤波器幅度分布);x=x1.*(h);subplot(3,1,3);plot(w/pi,abs(x);title(滤除噪声后的音频信号频谱图);结果分析要使噪声消除，则滤波器在噪声频率处需要完成衰减，此音频信号中，需要完成的衰减约为49.54db，当使用双线性变换法时，在此噪音频率前已已达到标准，过渡带衰减速率高。(4) 定量比较上述两种滤波器的各项性能:从第一题中可以看出，两方法的最大区别是过渡带的差距，则先保持ap=1;as=25不变，调大wp，使通带增加， np=6336w=linspace(0,2*pi

27、,50688);np=6336;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) ,blue); xla

28、bel(normalized frequency);ylabel(gain,db); hold on;np=6336;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;

29、plot(w/pi,20*log10(abs(h)/norm) ,red); xlabel(normalized frequency);ylabel(gain,db);legend(脉冲响应不变法, 双线性变化法)title(两种滤波器np=6336);np=6500w=linspace(0,2*pi,50688);np=6500;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);nu

30、ma,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , blue); xlabel(normalized frequency);ylabel(gain,db); hold on;np=6500;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(w

31、p,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脉冲响应不变法, 双线性变化法)title(两种滤波器np=6500);np=7000w=linspace(0,

32、2*pi,50688);np=7000;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , blue)

33、; xlabel(normalized frequency);ylabel(gain,db); hold on;np=7000;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/

34、norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脉冲响应不变法, 双线性变化法)title(两种滤波器np=7000);np=7500w=linspace(0,2*pi,50688);np=7500;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2

35、/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , blue); xlabel(normalized frequency);ylabel(gain,db); hold on;np=7200;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=but

36、tord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脉冲响应不变法, 双线性变化法)title(两种滤波器np=7200);np=7100w=linsp

37、ace(0,2*pi,50688);np=7100;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) ,

38、 blue); xlabel(normalized frequency);ylabel(gain,db); hold on;np=7100;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd

39、=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脉冲响应不变法, 双线性变化法)title(两种滤波器np=7100);np=8000w=linspace(0,2*pi,50688);np=8000;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-

40、1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=impinvar(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , blue); xlabel(normalized frequency);ylabel(gain,db); hold on;np=8000;ns=8448;n=50688;wp=np*2*pi/n; ws=ns*2*pi/n;ap=1;as=25;fs=1;wp=wp*fs; ws=ws*fs

41、;n=buttord(wp,ws,ap,as,s);wc=wp/(10(0.1*ap)-1)(1/2/n);numa,dena=butter(n,wc,s);numd,dend=bilinear(numa,dena,fs);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm) , red); xlabel(normalized frequency);ylabel(gain,db);legend(脉冲响应不变法, 双线性变化法)title(两种滤波器np=8000); 结果分析当ws，as,ap不变时，wp增大，两种方法的过渡带变小，阻带变宽，但是用脉冲响应法做出的变化剧烈。（5）hamming窗:（n=85）x,fs,bits=wavread(hello.wav); x1=x(1:1:end);x1=fft(x1); w=linspace(0,2*pi,50688);subplot(3,1,1);plot(w/pi,abs(x1);title(含有噪声的音频信号频谱图);np=6336;ns=8448;n=