IIR和FIR滤波器设计专题研讨.doc

数字信号处理课程研究性学习报告姓名 曾海宁 学号 11212023 同组成员 张欧(111212024) 易伟(11212022) 指导教师 时间 IIR和FIR滤波器设计专题研讨【目的】(1) 掌握根据滤波器指标设计IIR和FIR数字滤波器的原理和方法。(2) 熟悉通过IIR和FIR数字滤波器进行实际系统设计的方法。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】 1 设计一个数字滤波器，在频率低于的范围内，低通幅度特性为常数，并且不低于0.75dB。在频率和之间，阻带衰减至少为20dB。(1)试求满足这些条件的最低阶Butterworth滤波器。(2)试求满足这些条件的最低阶Chebyshev I滤波器。(3)自主选择一段带限信号，通过所设计的（1）、（2）两种滤波器，比较各自的输入和输出信号。讨论两种滤波器在结构和性能上的差异。【题目分析】本题讨论模拟滤波器和数字滤波器的设计。结合课本分析各类滤波器的特性。数字滤波器是一个离散系统，其系统函数一般可表示为z-1的有理多项式形式，即在设计数字滤波器时，首先根据具体的滤波成分，再确定设计数字滤波器的技术指标，再求出满足设计指标的离散系统的系统函数H(z)。题目中要求分别设计Butterworth滤波器和Chebyshev I滤波器，在选择一段带限信号通过这两个滤波器，实际上就要先将数字技术指标转换成模拟的技术指标，在根据该指标设计出模拟的低通滤波器，再转换成为所需的数字滤波器。BW的N CB1的N可以根据分式的关系来确定N的最小值【IIR模拟滤波器设计的基本方法】IIR滤波器设计的基本方法是首先将数字滤波器的设计指标转换成为模拟滤波器的设计指标,然后设计模拟滤波器H(s)时，再将模拟滤波器H(s)变换为数字滤波器H(z)。在将H(s)变换为H(z)时，要求模拟域到数字域的映射满足下列两个条件：（1） 两者的频率特性不变，即s平面的虚轴jw必须映射到z平面的单位圆上。（2） 变换后的滤波器仍然是稳定的，即s左半平面必须映射到z平面的单位圆内。这样才能保证变换后的数字滤波器的频率响应与模拟滤波器的频率响应基本一致。s到z域转换设计模拟滤波器频率转换待设计数字滤波器指标数字滤波器H(z)模拟滤波器H(s)模拟滤波器指标IIR数字滤波器的设计过程而模拟滤波器的设计都是通过设计模拟低通滤波器来实现，比较常用的模拟低通滤波器有Butterworth(巴特沃思)和Chebyshev(切比雪夫)等。将模拟滤波器变换为数字滤波器的主要方法有脉冲响应不变法和双线性变换法。【仿真结果】(1) Butterworth滤波器Ap= 0.4192 As= 19.9999 N = 8(2) Chebyshev I滤波器Ap= 0.7523 As= 21.3368 N = 4(3)1)信号通过 Butterworth滤波器前后的信号波形对比信号通过 Butterworth滤波器前后的频谱波形对比2) 信号通过 Chebyshev I滤波器前后的信号波形对比信号通过 Chebyshev I滤波器前后的频谱波形对比【结果分析】Butterworth滤波器通带拥有很好的稳定性，基本上恒定在一个常数上，然后整一个滤波器在频域范围内单调下降，衰减的也很快，是一个非常理想的滤波器。ChebysheyI滤波器在通带会有稍微的波动，而除去通带外的部分也是呈现单调下降的波形。【自主学习内容】（1）设计IIR数字滤波器的步骤。（2）用Butterworth、ChebyshevI和 ChebyshevII等设计方法设计模拟低通滤波器。（3）通过频率转换变换将模拟低通滤波器转换为所需的滤波器的类型。【阅读文献】数字信号处理（第二版） 陈后金主编 高等教育出版社【发现问题】 (专题研讨或相关知识点学习中发现的问题)：在自己设计的带限信号的高阶部分并不能被滤波器所滤除。【问题探究】【仿真程序】(1) 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);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);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= %.4f ,-20*log10( abs(h(1);fprintf(As= %.4f ,-20*log10( abs(h(2);disp(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,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 Ws;numd=numd/norm;h=freqz(numd,dend,w);fprintf(Ap= %.4f ,-20*log10( abs(h(1);fprintf(As= %.4f ,-20*log10( abs(h(2);disp(N = );disp(N);(3) 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);numa,dena=butter(N,wc,s);numd,dend=bilinear(numa,dena,fs);w=linspace(0,pi,1024);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;t=0:0.01:2*pi;y=sin(5*pi*t)+0.5*sin(50*pi*t);y1=filter(numd,dend,y);Y=fft(y);Y1=fft(y1);figure(1);subplot(2,1,1);plot(abs(Y);title(before(Butterworth);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(2,1,2);plot(abs(Y1);title(after(Butterworth);ylabel(Gain in dB);xlabel(Fs in Hz);figure(2);subplot(2,1,1);plot(y);title(before(Butterworth);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(2,1,2);plot(y1);title(after(Butterworth);ylabel(Gain in dB);xlabel(Fs in Hz);2)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);numa,dena=cheby1(N,Ap,wc,s);numd,dend=bilinear(numa,dena,fs);w=linspace(0,pi,1024);h=freqz(numd,dend,w);norm=max(abs(h);numd=numd/norm;t=0:0.01:2*pi;y=sin(5*pi*t)+0.5*sin(50*pi*t);y1=filter(numd,dend,y);Y=fft(y);Y1=fft(y1);figure(1);subplot(2,1,1);plot(abs(Y);title(before(Butterworth);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(2,1,2);plot(abs(Y1);title(after(Butterworth);ylabel(Gain in dB);xlabel(Fs in Hz);figure(2);subplot(2,1,1);plot(y);title(before(Butterworth);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(2,1,2);plot(y1);title(after(Butterworth);ylabel(Gain in dB);xlabel(Fs in Hz);2 分别用Hamming窗，Blackman窗和Kaiser窗设计，满足下列指标的FIR低通滤波器：(1)画出所设计滤波器的幅频响应。(2)比较这三种窗的设计结果。【题目分析】本题讨论窗函数法设计数字FIR滤波器。结合课本分析不同窗函数法的设计结果。窗函数法设计线性FIR数字滤波器时：矩形窗的主瓣宽度为4pi/N,Ap=0.82dB,As=21dB，所设计的FIR滤波器过渡带宽度近似为1.8pi/N； Hamming窗主瓣宽度为8pi/N，Ap=0.019dB,As=53dB,所设计的FIR滤波器过渡带宽度近似为6.2pi/N；Blackman窗主瓣宽度为12pi/N，Ap=0.0017dB,As=74dB,s所设计的FIR滤波器过渡带宽度近似为11.4pi/N；Kaiser窗可以通过调节 参数调整窗函数的形状，从而达到不同的阻带衰减。【FIR模拟滤波器设计的基本方法】窗函数法设计FIR滤波器的基本思想是在时域逼近理想滤波器的单位脉冲响应.首先根据待逼近的理想滤波器的频率响应,由IDTFT求出理想滤波器的单位脉冲响应,再将无限长的加窗截断得到有限长序列 。为了获得线性相位FIR滤波器，在窗函数法设计FIR滤波器哦的过程中，需要将线性相位因子 加入理想滤波器的频率响应。【仿真结果】1)Hamming窗N=352) Blackman窗N=593) Kaiser窗N=27【结果分析】（2）Hamming窗，Blackman窗和Kaiser窗三种窗中Kaiser窗的阶数最小，衰减也较快，Blackman窗衰减最大。可以看出，同一个窗函数的阶数越高，其衰减越大，但是，这样会增加计算量。Blackman窗的衰减最大，可是其阶数也是最高的。可见窗函数的阶数跟衰减速度也是一对矛盾。然后Kaiser窗可以在一定程度上改善这样的情况。【自主学习内容】(1) 窗函数设计线性相位FIR滤波器的基本步骤。(2) 使用MATLAB实现FIR滤波器。【阅读文献】数字信号处理（第二版） 陈后金主编 高等教育出版社【发现问题】 (专题研讨或相关知识点学习中发现的问题)：【问题探究】【仿真程序】(1)1)Hamming窗Wp=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);w=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);xlabel(hamming);grid on;2) Blackman窗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=%.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,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb);title(blackman);ylabel(Gain in dB);xlabel(F/Hz);grid on;3) Kaiser窗Ap=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);omega=linspace(0,pi,512);mag=freqz(h,1,omega);plot(omega/pi,20*log10(abs(mag);xlabel(Normalized frequency);title(Kaiser)ylabel(Gain in db);legend(Kaiser);grid on;3附件给出了一段含有噪声的音频信号。（1）分析该信号的频谱特点。（2）通过脉冲响应不变法设计一个滤波器对其进行处理，得到有用信息，自主确定各项指标。（3）使用双线性变换法重做（2）。（4）试定量比较上述两种滤波器的各项性能，画出能说明性能差异的相关图形，对比并解释。（5）通过不同窗函数法设计FIR滤波器对其进行处理，比较不同窗函数的滤波效果，给出分析和结论。（选做）(6) 请尝试采用其它的音频信号，混入不同的噪声，利用所学的滤波方法进行分析，会得到什么样的效果？【题目分析】本题讨论用IIR和FIR数字滤波器进行实际系统设计的方法。【仿真结果】(1)(2) 脉冲响应不变法Ap= 0.7959 As= 20.0000N = 17滤波器的幅频特性及相频特性滤波前后的频谱图滤波前后的信号的波形(3) 双线性变换法Ap= 3.4894As= 31.4938 N = 17滤波器的幅频特性及相频特性滤波前后的频谱图滤波前后的信号的波形(4)双线性变换法跟脉冲响应不变法的比较(5) 信号通过各个滤波器之后的频谱图输出信号之间的频谱比较【结果分析】对几种滤波器得到的音频信号进行对比。（1）信号的频谱总的是在低频，但是在中频部分有两个毛刺。并且幅值跟信号的幅值差不多，所以需要滤除。（2）设计低通滤波器进行对毛刺的滤除。脉冲响应不变法：Fp=1000;Fs=1200;Ft=8000; Ap=1; As=20；得到的Ap= 0.7959，As= 20.0000，N = 17。（3）设计低通滤波器进行对毛刺的滤除。双线性变换法：Fp=1000;Fs=1200;Ft=8000; Ap=1; As=20；得到的Ap= 3.4894，As= 31.4938，N = 17。（4）根据上面的仿真分析，可以看出，同样的滤波器设计要求，利用脉冲响应不变法设计出来的滤波器的滤波效果要比利用双线性变换法设计出来的滤波器哦的效果要差。但是，总体效果还是不错的。由上可以看出，通带截频的增益上看，脉冲响应不变法的效果要优于双线性变换法，但是阻带截频的增益上，则反过来了。（5）在上面的仿真分析，可以看出，在fp=100;fc=120;As=60;Ap=1;Fs=1000;的设计要求下，凯泽窗和blackman窗的滤波器效果相似，都比hamming窗的效果好得多。也就是说，在这样同等的低通滤波的要求下，凯泽窗和blackman窗的效果要好一些。（6）经过滤波之后的音频信号的噪声明显小了，但是仍然不及原音频信号的效果好。加入了增益为30db的高斯白噪声，经过滤波之后，噪声明显下降。【自主学习内容】（1）awgn的使用方法；（2）filter的使用方法；（3）fftfilt的使用方法；【阅读文献】数字信号处理（第二版） 陈后金主编 高等教育出版社【发现问题】 (专题研讨或相关知识点学习中发现的问题)：使用脉冲响应不变法设计滤波器时，由于阻带增益太大，导致输出出现错误。【问题探究】在IIR滤波器设计过程中，由于利用指标数据的不同，造成裕量的出现。讨论利用不同指标出现裕量对滤波器性能的影响，以及如何有效地利用它？答：由于裕量的存在，我们可以适当的降低滤波器的要求，充分利用裕量。例如当以通带增益为标准设计时，可以根据设计出来的滤波器的性能适当的减少阻带的增益。如果，裕量的大量存在，会导致滤波器的设计的不必要，增加计算机的计算量，降低计算速率。在FIR滤波器设计过程中，由于所选用窗函数的不同，导致对信号滤波的效果不同。本题的语音信号加入的是双频噪声，讨论如果对信号加入其它形式的噪声，采用何种滤波器更合适。答：【仿真程序】(1)y,fs,bits=wavread(Hello for discussion.wav);sound(y,fs,bits);N=length(y);Y=fft(y,N);subplot(2,1,2);plot(abs(Y);grid;title(Frequency domain);xlabel(Fs in Hz);subplot(2,1,1);plot(y);grid;xlabel(time);title(Time domain);(2) 脉冲响应不变法y,fs,bits=wavread(Hello for discussion);sound(y,fs,bits);Y=fft(y);Fp=1000;Fs=1200;Ft=8000;Wp=2*pi*Fp/Ft; Ws=2*pi*Fs/Ft; Ap=1; As=20;fsam=1; wp=Wp*fsam; ws=Ws*fsam;N,wc=buttord(wp,ws,Ap,As,s);numa,dena=butter(N,wc,s);numd,dend=impinvar(numa,dena,fsam);w=linspace(0,pi,1024);h=freqz(numd,dend,w);norm=max(abs(h);w=Wp Ws;numd=numd/norm;h=freqz(numd,dend,w);freqz(numd,dend);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2);y1=filter(numd,dend,y);Y1=fft(y1);figure(2);subplot(2,1,1);plot(abs(Y);title(before);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(2,1,2);plot(abs(Y1);title(after);ylabel(Gain in dB);xlabel(Fs in Hz);figure(3);subplot(2,1,1);plot(y);title(before);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(2,1,2);plot(y1);title(after);ylabel(Gain in dB);xlabel(Fs in Hz);sound(y1,fs,bits);(3) 双线性变换法y,fs,bits=wavread(Hello for discussion);sound(y,fs,bits);Y=fft(y);Fp=1000;Fs=1200;Ft=8000;Wp=2*pi*Fp/Ft; Ws=2*pi*Fs/Ft; Ap=1; As=20;fsam=1; wp=Wp*fsam; ws=Ws*fsam;N,wc=buttord(wp,ws,Ap,As,s);numa,dena=butter(N,wc,s);numd,dend=bilinear(numa,dena,fsam);w=linspace(0,pi,1024);h=freqz(numd,dend,w);norm=max(abs(h);w=Wp Ws;numd=numd/norm;h=freqz(numd,dend,w);freqz(numd,dend);fprintf(Ap= %.4fn,-20*log10( abs(h(1);fprintf(As= %.4fn,-20*log10( abs(h(2);y1=filter(numd,dend,y);Y1=fft(y1);figure(2);subplot(2,1,1);plot(abs(Y);title(before);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(2,1,2);plot(abs(Y1);title(after);ylabel(Gain in dB);xlabel(Fs in Hz);figure(3);subplot(2,1,1);plot(y);title(before);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(2,1,2);plot(y1);title(after);ylabel(Gain in dB);xlabel(Fs in Hz);sound(y1,fs,bits);(4)y,fs,bits=wavread(Hello for discussion);Y=fft(y);%Fp=1000;Fs=1200;Ft=8000;Wp=2*pi*Fp/Ft; Ws=2*pi*Fs/Ft; Ap=1; As=20;fsam=1; wp=Wp*fsam; ws=Ws*fsam;N,wc=buttord(wp,ws,Ap,As,s);num,den=butter(N,wc,s);numd,dend=impinvar(num,den,fsam);%numa,dena=bilinear(num,den,fsam);%w=linspace(0,pi,1024);h=freqz(numd,dend,w);h2=freqz(numa,dena,w);norm=max(abs(h);norm2=max(abs(h2);w=Wp Ws;numd=numd/norm;numa=numa/norm2;y1=filter(numd,dend,y);y2=filter(numa,dena,y);%Y1=fft(y1);Y2=fft(y2);figure(1);subplot(2,1,1);plot(abs(Y1)-abs(Y2);title(impinvar - bilinear);xlabel(Fs in Hz);ylabel(Gain in dB);subplot(2,1,2);plot(abs(Y);title(orginal);xlabel(Fs in Hz);ylabel(Gain in dB);(5) y,fs,bits=wavread(Hello for discussion);Y=fft(y);fp=100;fc=120;As=60;Ap=1;Fs=1000;wc=2*pi*fc/Fs; wp=2*pi*fp/Fs;wd=abs(wc-wp);beta=0.1102*(As-8.7);%?betaN=ceil(As-7.95)/2.285/wd);%?Nwn= kaiser(N+1,beta);ws=(wp+wc)/2/pi;Kwin=fir1(N,ws,wn);%?N=ceil(7*pi/wd);N=mod(N+1,2)+N;%HammingNwn= hamming(N);Hwin=fir1(N-1,ws,wn);%HammingN=ceil(11.4*pi/wd);N=mod(N+1,2)+N;%BlackmanNwn= blackman(N);Bwin=fir1(N-1,ws,wn);%Blackmanb=fftfilt(Bwin,y);B=fft(b);h=fftfilt(Hwin,y);H=fft(h);k=fftfilt(Kwin,y);K=fft(k);figure(1);subplot(3,1,1);plot(abs(B)-abs(H);grid on;title(blackman - hamming);ylabel(Gain in dB);xlabel(Fs in Hz);subplot(3,1,2);plot(abs(K)-abs(H);grid on;title(kaiser - hamming);ylabel(Gain in dB);xlabel(Fs in Hz)