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

1、数字信号处理课程研究性学习报告IIR和FIR滤波器设计专题研讨【目的】(1) 掌握根据滤波器指标设计IIR和FIR数字滤波器的原理和方法。(2) 熟悉通过IIR和FIR数字滤波器进行实际系统设计的方法。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】 1 设计一个数字滤波器，在频率低于的范围内，低通幅度特性近似为常数，并且不低于0.75dB。在频率和之间，阻带衰减至少为20dB。(1)试求满足这些条件的最低阶Butterworth滤波器。(2)试求满足这些条件的最低阶Chebyshev I滤波器。(3)自主选择一段带限信号，通过所设计的（1）、（2）两种滤波器

2、，比较各自的输入和输出信号。讨论两种滤波器在结构和性能上的差异。【题目分析】本题讨论模拟滤波器和数字滤波器的设计。结合课本分析各类滤波器的特性。【IIR模拟滤波器设计的基本方法】【仿真结果】1)N=8N=2N=12)N=4N=2N=1【结果分析】BW的N CB1的N可以根据分式的关系来确定N的最小值BW的N和Ap成正比CB1的N和Wp成反比在对其取值的时候可以根据这样的规律实验，模拟的结果和理论是相同的。从图像上可以看出N的值就是通频带的起伏的次数，N=1就起伏一次，N=2就起伏2次，N=8就接近连成一条线。【自主学习内容】【阅读文献】【仿真程序】1)Wp=0.2613*pi; Ws=0.40

3、18*pi; Ap=0.75; As=20;Fs=1; wp=Wp*Fs; ws=Ws*Fs;N=buttord(wp,ws,Ap,As,s);wc=wp/(10(0.1*Ap)-1)(1/N/2);numa,dena=butter(N,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h/norm);xlabel( butter );ylabel(Ap=0.75,dB);d

4、isp(N=);disp(N);2)Wp=0.2613*pi; Ws=0.4018*pi; Ap=0.75; As=20;Fs=1; wp=Wp*Fs; ws=Ws*Fs;N=cheb1ord(wp,ws,Ap,As,s);wc=wp/(10(0.1*Ap)-1)(1/N/2);numa,dena=cheby1(N,Ap,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h/

5、norm);xlabel(cheby1);ylabel(Ap=0.75,dB);disp(N=);disp(N);2 分别用Hamming窗，Blackman窗和Kaiser窗设计，满足下列指标的FIR低通滤波器：(1)画出所设计滤波器的幅频响应。(2)比较这三种窗的设计结果。【题目分析】本题讨论窗函数法设计数字FIR滤波器。结合课本分析不同窗函数法的设计结果。【FIR模拟滤波器设计的基本方法】【仿真结果】【结果分析】【仿真程序】Wp=0.4*pi;Ws=0.6*pi;Ap=0.5;As=45;%hanningfigure(1)N=floor(6.2*pi/(Ws-Wp);N=mod(N+1,

6、2)+N;M=N-1;w=hanning(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);plot(omega/pi,20*log10(abs(mag);grid onxlabel(Normalized frequency)ylabel(Gain in dB)title(hanning) %hammingfigure(2)N=floor(7*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1w=hamming(N

7、);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);plot(omega/pi,20*log10(abs(mag);grid onxlabel(Normalized frequency)ylabel(Gain in dB)title(hamming) %blackmanfigure(3)N=floor(11.4*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1;w=blackman(N);Wc=(Wp+Ws)/2;k

8、=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);plot(omega/pi,20*log10(abs(mag);grid onxlabel(Normalized frequency)ylabel(Gain in dB)title(blackman) %kaiserfigure(4)f=Wp/pi Ws/pi;a=1,0;dev=10(-As/20)*ones(1,length(a);M,Wc,beta,ftype=kaiserord(f,a,dev);M=m

9、od(M,2)+M;h=fir1(M,Wc,ftype,kaiser(M+1,beta);omega=linspace(0,pi,512);mag=freqz(h,1,omega);plot(omega/pi,20*log10(abs(mag);grid onzoom onxlabel(Normalized frequency)ylabel(Gain in dB)title(kaiser)3附件给出了一段含有噪声的音频信号。（1）分析该信号的频谱特点。（2）通过脉冲响应不变法设计一个滤波器对其进行处理，得到有用信息，自主确定各项指标。（3）使用双线性变换法重做（2）。（4）试定量比较上述两种滤

10、波器的各项性能，画出能说明性能差异的相关图形，对比并解释。（5）通过不同窗函数法设计FIR滤波器对其进行处理，比较不同窗函数的滤波效果，给出分析和结论。（选做）(6) 请尝试采用其它的音频信号，混入不同的噪声，利用所学的滤波方法进行分析，会得到什么样的效果？【题目分析】本题讨论用IIR和FIR数字滤波器进行实际系统设计的方法。【仿真结果】z1= wavread(C:UserszhangDesktop课件DSP你好加噪声信号);plot(z1);z1= wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);n=0:81

11、91;plot(n,Y1);（1）窗函数设计低通滤波器 程序设计如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1;Fs=8000;wc=2*pi*fc/Fs; wp=2*pi*fp/Fs;wdel=wc-wp;beta=0.112*(As-8.7);N=ceil(As-8)/2.285/wdel);wn= kaiser(N+1,beta); ws=(wp+wc)/2/pi;b=fir1(N,ws,wn

12、);figure(1);freqz(b,1);x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波后信号波形);sound(x,fs,bits);wavplay(x,44100);wav

13、write(x,44100,1b.wav);fname=sprintf(E:1b.wav); wavwrite(x,44100,fname);图形分析如下：（2）窗函数设计高通滤波器程序设计如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);fp=2800;fc=3000;As=100;Ap=1;Fs=8000;wc=2*pi*fc/Fs; wp=2*pi*fp/Fs;wdel=wc-wp;beta=0.112*(As-8.7);N=ceil(As-8)/2

14、.285/wdel);wn= kaiser(N,beta); ws=(wp+wc)/2/pi;b=fir1(N-1,ws,high,wn);figure(1);freqz(b,1);x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot

15、(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,2b.wav);fname=sprintf(E:2b.wav); wavwrite(x,44100,fname);（3）窗函数设计带通滤波器程序设计如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);fp1=1200 ;fp2=3000 ;fc1=1000 ;fc2=3200 ;As=100 ;A

16、p=1 ;Fs=8000 ;wp1=2*pi*fp1/Fs; wc1=2*pi*fc1/Fs; wp2=2*pi*fp2/Fs; wc2=2*pi*fc2/Fs;wdel=wp1-wc1;beta=0.112*(As-8.7);N=ceil(As-8)/2.285/wdel);ws =(wp1+wc1)/2/pi,(wp2+wc2)/2/pi;wn= kaiser(N+1,beta); b=fir1(N,ws,wn);figure(1);freqz(b,1)x=fftfilt(b,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);a

17、xis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,2000,0,0.0003);title(滤波后信号频谱)subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,3b.wav);fname=sprintf(E:3b.wav); wavwrite(x,44100,fname);图形分析如下：（4）双线性变换法设计低

18、通滤波器选用butter程序设计如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1;Fs=8000;wc=2*fc/Fs; wp=2*fp/Fs; N,ws=buttord(wc,wp,Ap,As);b,a=butter(N,ws);figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);pl

19、ot(abs(Y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,4b.wav);fname=sprintf(E:4b.wav); wavwrite(x,44100,fname);图形分析选用che

20、by1程序设计如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);fp=1000;fc=1200;As=100;Ap=1; ;Fs=8000;wc=2*fc/Fs;wb=2*fp/Fs;n,wp=cheb1ord(wc,wb,Ap,As);b,a=cheby1(n,Ap,wp);figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);a

21、xis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,4000,0,0.03);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,5b.wav);fname=sprintf(E:5b.wav); wavwrite(x,44100,fname);图形分析如下：5），双线性变换法设计高通

22、滤波器选用butter程序设计如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);fc=2800 ;fp=3000 ;As=100;Ap=1;Fs=8000;wc=2*fc/Fs; wp=2*fp/Fs; N,ws=buttord(wc,wp,Ap,As);b,a=butter(N,ws,high);figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2

23、,1);plot(abs(Y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,6b.wav);fname=sprintf(E:6b.wav); wavwrite(x,44100,fname);图形分

24、析如下：选用cheby1程序设计如下：clear;close allz1,fs,bits=wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);fc=2800 ;fp=3000 ;As=100;Ap=1; Fs=8000;wc=2*fc/Fs;wb=2*fp/Fs;n,wp=cheb1ord(wc,wb,Ap,As);b,a=cheby1(n,Ap,wp,high);figure(1);freqz(b,a);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1)

25、;plot(abs(Y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,4000,0,0.03);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);wavplay(x,44100);wavwrite(x,44100,7b.wav);fname=sprintf(E:7b.wav); wavwrite(x,44100,fname);图形分析如

26、下：（6），双线性变换法设计带通滤波器选用butter程序设计如下clear;close allz1,fs,bits=wavread(C:UserszhangDesktop课件DSP你好加噪声信号);y1=z1(1:8192);Y1=fft(y1);fp1=1200 ;fp2=3000; fc1=1000 ;fc2=3200 ;As=100;Ap=1; Fs=8000; wc=2*fc1/Fs,2* fc2/Fs;wp=2*fp1/Fs,2*fp2/Fs; N,ws=buttord(wc,wp,Ap,As);b,a=butter(N,ws,stop);figure(1);freqz(b,a,512,Fs);x=filter(b,a,z1);X=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(Y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,100