FIR和IIR滤波器设计

1、IIR和FIR滤波器设计专题研讨第八组：彭然 09214040 海骁飞09214032 季凯09214035 刘炘09214036 1 设计一个数字滤波器，在频率低于 的范围内，低通幅度特性近似为常数，并且不低于0.75dB。在频率 和 之间，阻带衰减至少为20dB。(1)试求满足这些条件的最低阶Butterworth滤波器。(2)试求满足这些条件的最低阶Chebyshev I滤波器。(3)自主选择一段带限信号，通过所设计的（1）、（2）两种滤波器，比较各自的输入和输出信号。讨论两种滤波器在结构和性能上的差异。0.2613 0.4018 【题目分析题目分析】本题讨论模拟滤波器和数字滤波器的设计

2、。结合课本分析各类滤波器的特性。分析：Butterworth滤波器幅度响应在通带和阻带都是随频率单调变化，在通带和阻带都存在裕量，并且阶数较高；Chebyshev I滤波器在通带具有等波纹特性，在阻带单调下降，因此通带裕量较Butterworth滤波器要小，但是阻带较Butterworth滤波器要大，阶数要低。【IIR模拟滤波器设计的基本方法模拟滤波器设计的基本方法】模拟滤波器的设计就是求出满足技术指标的连续系统函数H（s）,模拟滤波器的技术指标有通带截频wp、通带最大衰减Ap、阻带截频ws和阻带最小衰减As。在设计模拟滤波器时，先将待设计的模拟滤波器技术指标转换为模拟低通滤波器技术指标，然后

3、设计模拟低通滤波器，再通过频率变换将模拟低通滤波器转换为所需的滤波器。其中模拟低通滤波器的设计有一些模型，比如Butterworth、Chebyshev等等模拟滤波器模型。（1）Butterworth滤波器:%ButterworthWp=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);numa,dena=butter(N,wc,s);numd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,1024);h=freqz(

4、numd,dend,w);norm=max(abs(h);plot(w/pi,20*log10(abs(h)/norm);xlabel(Normalized frequency);ylabel(Gain,dB);grid on;w=Wp Ws;numd=numd/norm;h=freqz(numd,dend,w);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2); N=8 Ap= 0.4192 As= 19.9999(2) Chebyshev I滤波器:%Chebyshev IWp=0.261

5、3*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);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);plot(w/pi,20*log10(abs(h)/norm);xlabel(Normalized frequency);ylabel(Gain,dB);grid on;w=Wp W

6、s;numd=numd/norm;h=freqz(numd,dend,w);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2);N=4 Ap= 0.7523 As= 21.3368（3）将带限信号x=cos(0.5*pi*t)通过Butterworth滤波器观察滤波效果：%ButterworthWp=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);numa,d

7、ena=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;H=freqz(numd,dend,w);t=0:0.01:2*pi;x=cos(0.5*pi*t);x1=filter(numd,dend,x);subplot(2,1,1);plot(t,x);hold on;plot(t,x1,r);title(Butterworth);grid on;subplot(2,1,2);plot(t,x1-x)

8、;title();grid on;将带限信号x=cos(0.5*pi*t)通过Chebyshev I滤波器观察滤波效果：%Chebyshev IWp=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);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/n

9、orm;H=freqz(numd,dend,w);t=0:0.01:2*pi;x=cos(0.5*pi*t);x1=filter(numd,dend,x);subplot(2,1,1);plot(t,x);hold on;plot(t,x1,r);title(Chebyshev I);grid on;subplot(2,1,2);plot(t,x1-x);title();grid on;【结果分析结果分析】(1) Butterworth滤波器幅度响应在通带和阻带都是随频率单调变化，在通带和阻带都存在裕量，设计要求是Ap=0.75，As=20，而结果是Ap= 0.4192 As= 19.9999

10、，并且阶数较高，N=4；（2）Chebyshev I滤波器在通带具有等波纹特性，在阻带单调下降，因此通带裕量较Butterworth滤波器要小，但是阻带较Butterworth滤波器要大，但是结果是Ap= 0.7495 As= 25.5334，通带衰减不满足题目要求，阻带存在较大的裕量。（3）通过分析滤波器后时域误差可知：Butterworth滤波器产生的误差大概是0.08498，要小于Chebyshev I滤波器的0.096。同时可知，时域波形在0处误差几乎是1，我觉得这个是由于余弦信号时域直接加了矩形窗，在频域为二者卷积，可看成是两个sa函数形式的波，其幅值与矩形窗长度有关，其N越大，幅值

11、越大，那么即使在滤波时幅值下降了，那么这个可以补偿，其时域也几乎一样。当t=0:0.01:2*pi;x=cos(0.5*pi*t);误差上图可知，当t=0:0.001:2*pi;x=cos(0.5*pi*t);2 分别用Hamming窗，Blackman窗和Kaiser窗设计，满足下列指标的FIR低通滤波器：0.4 ,0.6 ,0.5,45pspsAdBAdB (1)画出所设计滤波器的幅频响应。(2)比较这三种窗的设计结果。【FIR模拟滤波器设计的基本方法模拟滤波器设计的基本方法】FIR数字滤波器的设计通常是根据理想滤波器的频率响应 ，采用窗函数法、频率取样法或优化方法使设计的滤波器的频率响应

12、 逼近 。jH ejH ejdHe%hammingFIRWp=0.4*pi;Ws=0.6*pi;Ap=0.5;As=45;N=ceil(7*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1;fprintf(N=%.0fn,N);%generate windoww=hamming(N);Wc=(Wp+Ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M);h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb)

13、;xlabel(hammingFIR);grid on; （1）Hamming窗：（1）Hamming窗：N=35 As=55（2）Blackman窗：%blackmanFIRWp=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=%.0fn,N);w=blackman(N);Wc=(Wp+Ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M);h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,ome

14、ga);magdb=20*log10(abs(mag);plot(omega/pi,magdb);xlabel(blackmanFIR);grid on;（2）Blackman窗：N=27 As=45%Program:KaiserFIRAp=0.5;As=45;Rp=1-10.(-0.05*Ap);Rs=10.(-0.05*As);f=0.4,0.6;a=1,0;dev=Rp,Rs;M,Wc,beta,ftype = kaiserord(f,a,dev);M=mod(M,2)+M;fprintf(N=%.0fn,M+1);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);xlabel(KaiserFIR);grid on;（3）Kaiser窗：（3）Kaiser窗：N=59 As=78.11【结果分析结果分析】Hamming窗，Blackman窗和Kaiser窗三种窗中Blackman窗的阶数最小，衰减也较快，Kaiser窗衰减最大3附件给出了一段含有噪声的音频信号。（1）分析该信号的频谱特点。（2）通过脉冲响应不变法设计一个滤波器对其进行处理，得到有用信息，自主确定各项指标。（3）使用双线性变换法重做（

16、2）。（4）试定量比较上述两种滤波器的各项性能，画出能说明性能差异的相关图形，对比并解释。（5）通过不同窗函数法设计FI R滤波器对其进行处理，比较不同窗函数的滤波效果，给出分析和结论。（选做）(6) 请尝试采用其它的音频信号，混入不同的噪声，利用所学的滤波方法进行分析，会得到什么样的效果？y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);%subplot(2,1,2);plot(abs(Y);xlabel();subplot(2,1,1);plot(y);xlabel();（1）信号的时域波形和频谱通过对频谱分析确定设计低通滤波器，其中通带截止

17、频率Fp=1400，阻带截止频率是Fs=1600Ap=1; As=20，用脉冲响应不变法：y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;Wp=2*pi*Fp/Ft; Ws=2*pi*Fs/Ft; Ap=1; As=20;fs=1; wp=Wp*fs; ws=Ws*fs;N,wc=buttord(wp,ws,Ap,As,s);numa,dena=butter(N,wc,s);numd,dend=impinvar(numa,dena,fs);w=linspace(0,pi,1024);h=f

18、reqz(numd,dend,w);norm=max(abs(h);w=Wp Ws;numd=numd/norm;h=freqz(numd,dend,w);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2);y1=filter(numd,dend,y);figure(1)subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);title(

19、);grid;subplot(2,2,1);plot(abs(Y),r);title();grid;（3）双线性变化法：y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Wp=2*pi*Fp/Ft;Ws=2*pi*Fs/Ft;wp=2*Ft*tan(Wp/2);ws=2*Ft*tan(Ws/2);N,wc=buttord(wp,ws,Ap,As,s); %num,den=butter(N,wc,s); %S numd,dend=bilinear(

20、num,den,Ft); %SZ h,w=freqz(numd,dend); %w=Wp Ws;h=freqz(numd,dend,w);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2);y1=filter(numd,dend,y);figure(1)subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);title();grid;su

21、bplot(2,2,1);plot(abs(Y),r);title();grid;（4）脉冲响应不变法与双线性变换法性能比较y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;Wp1=2*pi*Fp/Ft; Ws1=2*pi*Fs/Ft; Ap=1; As=20;fs=1; wp1=Wp1*fs; ws1=Ws1*fs;N1,wc1=buttord(wp1,ws1,Ap,As,s);numa1,dena1=butter(N1,wc1,s);numd1,dend1=impinvar(numa1,

22、dena1,fs);w=linspace(0,pi,1024);h1=freqz(numd1,dend1,w);norm=max(abs(h1);numd1=numd1/norm;h1=freqz(numd1,dend1,w);Wp2=2*pi*Fp/Ft;Ws2=2*pi*Fs/Ft;wp2=2*Ft*tan(Wp2/2);ws2=2*Ft*tan(Ws2/2);N2,wc2=buttord(wp2,ws2,Ap,As,s); num2,den2=butter(N2,wc2,s); numd2,dend2=bilinear(num2,den2,Ft); h2,w=freqz(numd2,de

23、nd2); y1=filter(numd1,dend1,y);y2=filter(numd2,dend2,y);plot(abs(y1-y2);xlabel();（5）通过不同窗函数法设计FI R滤波器处理：Hamming窗%HammingFIRy,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Wp=2*pi*Fp/Ft;Ws=2*pi*Fs/Ft;N=ceil(7*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1;fprintf(

24、N=%.0fn,N);%generate windoww=hamming(N);Wc=(Wp+Ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M);h=hd.*w;y1=fftfilt(h,y);figure(1)subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);title();grid;subplot(2,2,1);plot(abs(Y),r);title();grid;Blackman窗%B

25、lackmanFIRy,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Wp=2*pi*Fp/Ft;Ws=2*pi*Fs/Ft;N=ceil(11.4*pi/(Ws-Wp);N=mod(N+1,2)+N;M=N-1;fprintf(N=%.0fn,N);w=blackman(N);Wc=(Wp+Ws)/2;k=0:M;hd=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M);h=hd.*w;y1=fftfilt(h,y);figure(1)

26、subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);title();grid;subplot(2,2,1);plot(abs(Y),r);title();grid;Kaiser窗%KaiserFIRy,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Rp=1-10.(-0.05*Ap);

27、Rs=10.(-0.05*As);f=2800/16000,3200/16000;a=1,0;dev=Rp,Rs;M,Wc,beta,ftype = kaiserord(f,a,dev);M=mod(M,2)+M;fprintf(N=%.0fn,M+1);h = fir1(M,Wc,ftype,kaiser(M+1,beta);y1=fftfilt(h,y);figure(1)subplot(2,2,3);plot(abs(y1),r);title();grid;subplot(2,2,4);plot(y1,b);title();grid;subplot(2,2,2);plot(y,b);ti

28、tle();grid;subplot(2,2,1);plot(abs(Y),r);title();grid;性能对比：y,fs,bits=wavread(.wav);%N=length(y);%Y=fft(y,N);% Fp=1400;%Fs=1600;%Ft=8000;%As=20;%As=20Ap=1;%Ap=1Wp=2*pi*Fp/Ft;Ws=2*pi*Fs/Ft;N1=ceil(7*pi/(Ws-Wp);N1=mod(N1+1,2)+N1;M1=N1-1;fprintf(N1=%.0fn,N1);%generate windoww1=hamming(N1);Wc=(Wp+Ws)/2;k

29、=0:M1;hd1=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M1);h1=hd1.*w1;Rp=1-10.(-0.05*Ap);Rs=10.(-0.05*As);f=2800/16000,3200/16000;a=1,0;dev=Rp,Rs;M2,Wc,beta,ftype = kaiserord(f,a,dev);M2=mod(M2,2)+M2;fprintf(N2=%.0fn,M2+1);h2= fir1(M2,Wc,ftype,kaiser(M2+1,beta); N3=ceil(11.4*pi/(Ws-Wp);N3=mod(N3+1,2)+N3;M3=N3-1;fprintf(N3=%.0fn,N3);w3=blackman(N3);Wc=(Wp+Ws)/2;k=0:M3;hd3=(Wc/pi)*sinc(Wc/pi)*(k-0.5*M3);h3=hd3.*w3; y1=fftfilt(h1,y);