第4章___数字滤波器设计

1、第四章 数字滤波 第二节 数字滤波的设计通信与电子信息当中，在对信号作分析与处理时，常会用到有用信号叠加无用噪声的问题。这些噪声信号有的是与信号同时产生的，有的是在传输过程中混入的，在接收的信号中，必须消除或减弱噪声干扰，这是信号处理中十分重要的问题。根据有用信号与噪声的不同特性，消除或减弱噪声，提取有用信号的过程就称为滤波。滤波器的种类很多，实现方法也多种多样，本章利用Matlab来进行数字滤波器的设计。数字滤波器是一离散时间系统，它对输入序列x(n)进行加工处理后，输出序列y(n)，并使y(n)的频谱与x(n)的频谱相比发生某种变化。由DSP理论得知，无限长冲激响应（IIR）需要递归模型来

2、实现，有限长冲激响应（FIR）滤波器可以采用递归的方式也可采用非递归的方式实现。本章把FIR与IIR滤波器分别用Matlab进行分析与设计。数字滤波器的结构参看数字信号处理一书。数字滤波器的设计一般经过三个步骤：1给出所需滤波器的技术指标。2设计一个H（Z），使其逼近所需要的技术指标。3实现所设计的H（Z）。4.2.1 IIR数字滤波器设计设计IIR数字滤波器的任务就是寻求一个因果、物理可实现的系统函数H(z)，使它的频响H(ej)满足所希望得到的低通频域指标，即通带衰减Ap、阻带衰减Ar、通带截频p、阻带截频r。而其它形式的滤波器由低通的变化得到。 采用间接法设计IIR数字滤波器就是按给定的

3、指标，先设计一个模拟滤波器，进而通过模拟域与数字域的变换，求得物理可实现的数字滤波器。从模拟滤波器变换到数字滤波器常用的有：脉冲响应不变法和双线性变换法。IIR滤波器的设计过程如下数字频域指标模拟频域指标设计模拟滤波器H(S) 设计数字滤波器H(z)1. 模拟滤波器简介 模拟滤波器的设计方法已经发展得十分成熟，常用的高性能模拟低通滤波器有巴特沃斯型、切比雪夫型和椭圆型，而高通、带通、带阻滤波器则可以通过对低通进行频率变换来求得。必须指出，这三种滤波器都是非线性的相频特性。具体设计方法参见数字信号处理一书。 （1）巴特沃斯低通滤波器的特点是：通、阻带均为单调下降。这种单调下降的特性使得系统的误差

4、分布不均匀。在设计中，如果在通带满足指标，阻带指标就过于好了。或者说，在阻带满足指标，通带指标也过于好了。这样不利于以最小的阶数来满足设计指标。 （2）切比雪夫模拟低通滤波器 分成型和型：型是通带等波纹，阻带单调下降； 型滤波器是通带单调下降，阻带等波纹。切比雪夫的等波纹特性使得可以用较小的阶数设计出符合指标的滤波器。 （3）椭圆滤波器特点是通、阻带均为等波纹。椭圆滤波器的误差均匀分布在通、阻带，比较上面几种滤波器，在同样衰耗指标下，椭圆滤波器所要求的阶数最小。在同样阶数情况下，椭圆滤波器的通带到阻带的变化最陡峭、性能最好。 2模拟滤波器映射到数字滤波器 从模拟滤波器变换到数字滤波器常用的有：

5、冲激响应不变法和双线性变换法。 （1）冲激响应不变法冲激响应不变法的基本思路为：设模拟滤波器的单位冲激响应为ha（t），数字滤波器的单位脉冲响应h(n)，是对ha（t）的均匀取样，即 H(S)h(t) h(n)= h(nT) H(z) 冲激响应不变法由于保持了模拟滤波器的冲激响应在取样时刻的取样值，所以具有较好的时间特性。但由于从模拟域（s域）变换到数字域（z域）的映射关系不是一一对应的关系，因而频谱混叠难以避免, 所以在实际应用中有它的局限性。3 Matlab滤波器设计函数简介 使用MATLAB信号处理工具箱中提供的函数可以很容易地实现IIR滤波器设计。常用的IIR数字滤波器设计函数有： （

6、1）IIR滤波器阶数（order）选择buttord - 巴特沃斯（Butterworth）滤波器阶数选择 cheb1ord - 切比雪夫（Chebyshev）I型滤波器阶数选择 cheb2ord - 切比雪夫（Chebyshev）II型滤波器阶数选择 ellipord - 椭圆（Elliptic）滤波器阶数选择 buttord用来确定巴特沃斯滤波器的阶数： N,Wc = buttord (Wp,Wr,Ap,Ar,options) Ap,Ar为通带最大衰减和阻带最小衰减，以dB为单位； Wp,Wr为归一化通带截频和阻带截频。 Options=low, high, bandpass, stop,

7、 s, 分别对应低通、高通、带通、带阻、模拟滤波器，默认情况下，为低通或带通。对于带通和带阻滤波器Wp=Wp1,Wp2，Wp1和 Wp2为通带的上、下截频，Wr=Wr1,Wr2，Wr1和Wr2分别为阻带的上、下截频。若为模拟滤波器定阶，令Options= s, 则Wp和Wr以（归一或非归一）给出，因此可大于1。输出参数N为滤波器的阶数，Wc为归一化3dB截频。若Ap=3dB，则Wc= Wp。 其他滤波器的定阶函数类似buttord的格式： 切比雪夫I型滤波器N,Wc = cheb1ord(Wp, Wr, Ap, Ar) 切比雪夫II型滤波器N,Wc = cheb2ord(Wp, Wr, Ap,

8、 Ar) 椭圆滤波器N,Wc = ellipord(Wp, Wr, Ap, Ar) （2）完整的滤波器设计函数，根据滤波器阶数和通带截频设计出实际数字滤波器的H(z) （模拟域到数字域的变换只采用双线性变换法）： butter - 巴特沃斯（Butterworth）滤波器设计 cheby1 -切比雪夫（Chebyshev）I型滤波器设计 cheby2 - -切比雪夫（Chebyshev）II型滤波器设计 ellip - 椭圆（Elliptic）滤波器设计 maxflat - 通用的巴特沃斯（Butterworth）低通滤波器设计 yulewalk - Yule-Walker滤波器设计（直接数字

9、滤波器设计法） butter的格式如下： b,a = butter(N,Wc,options) z,p,k = butter(N,Wc,options) A,B,C,D = butter(N,Wc,options) N为滤波器的阶数，Wc为频率归一化后的3dB截频：0 < Wc < 1。 Options=low, high, bandpass, stop, s, 分别对应低通、高通、带通、带阻、模拟滤波器，默认情况下，为低通或带通。对于带通和带阻滤波器，Wc=W1,W2，W1和W2分别为通带的上、下截频。 输出变量的三种格式对应着滤波器系统函数H (z)的分子分母多项式形式、零极点

10、-增益形式和状态空间形式。若设计模拟滤波器，令Options= s, 则Wc以rad/ s（归一化或非归一化）给出，因此可大于1。例1 设计一个数字低通滤波器，要求在00.2内衰耗不大于3dB，在0.6内衰耗不小于60dB。 解：若模拟低通原型采用巴特沃斯滤波器： N,Wc=buttord(0.2,0.6,3,60); b,a=butter(N,Wc); freqz(b,a);axis(0,1,-120,0); 可见，相位为非线性。 例2 使用脉冲响应不变法设计一个数字带通滤器，其指标为：1）3dB通带中心频率在 0.52）上截频 = 0.55,下截频 0.453）在阻带 0.6处，最小衰耗A

11、10dB设取样周期T5×106s,要求>衰耗曲线单调增加。解：（1）确定数字频率指标 omega0=0.5*pi; %通带中心频率 omegac1=0.45*pi ; omegac2=0.55*pi; Apmax=3; %通带截频 omegar2=0.6*pi; Armin=10; T=5*pi*10(-6); %取样间隔（2）利用脉冲响应不变法设计数字滤波器1）按给定数字频率求得相应的模拟角频率wc1=omegac1/T; wc2=omegac2/T;w0=omega0/T;omegar1=omega0-(omegar2-omega0);wr1=omegar1/T; wr2=

12、omegar2/T;B=wc2-wc1;2)将指标转换为归一化模拟低通频率指标norm_wr1=(wr12)-(w02)/(B*wr1);norm_wr2=(wr22)-(w02)/(B*wr2);norm_wc1=(wc12)-(w02)/(B*wc1);norm_wc2=(wc22)-(w02)/(B*wc2);3)设计归一化模拟低通滤波器确定归一化截频if (abs(norm_wr1)-abs(norm_wr2)>0 norm_wr=abs(norm_wr2);else norm_wr=abs(norm_wr1);endnorm_wc=1;确定归一化模拟低通滤波器阶数 N=butt

13、ord(norm_wc,norm_wr,Apmax ,Armin,'s');设计归一化模拟低通滤波器 b_LP,a_LP=butter(N,norm_wc,'s');4)把归一化模拟低通滤波器转换为模拟带通滤波器b_BP,a_BP=lp2bp(b_LP,a_LP,w0, B);5)利用脉冲响应不变法把模拟滤波器转换为数字滤波器Fs=1/T;b,a=impinvar(b_BP,a_BP,Fs);w=linspace(0,2*pi,500);h=freqz(b,a,w);plot(w,20*log10(abs(h); axis(0,2*pi,-50,20); gri

14、dxlabel('frequency (rad)');ylabel('gain (dB)');sys=tf(b,a,T)(2) 双线性变换法 双线性变换法克服了脉冲响应不变法存在频谱混叠的缺点, 但在设计过程中一定要注意，由于宽频带压缩到较窄的频率范围，导致了频率畸变。其基本思想是：首先按给定的指标设计一个模拟滤波器，通过适当的数学变换方法将s域映射到z域，把无限宽的频带变换成有限宽频带。由于在实现数字化以前已经对频带进行了压缩，所以在数字化以后的频响可以做到无频谱混叠。该法可以用来设计各种不同频带要求的数字滤波器，所以获得广泛应用。例3：按例2的指标使用双线性

15、变换法设计一个数字带通滤波器。解： （1）给定数字频率指标omega0=0.5*pi; %通带中心频率omegac1=0.45*pi; omegac2=0.55*pi; Apmax=3; %通带截频omegar2=0.6*pi; Armin=10;T=5*pi*10(-6); %取样间隔（2）用双线性变换法设计滤波器1）按给定数字频率求得相应的模拟角频率：频率预畸变omegar1=omega0-(omegar2 - omega0);wc1=(2/T)*tan(omegac1/2); wc2=(2/T)*tan(omegac2/2);wr1=(2/T)*tan(omegar1/2); wr2=(

16、2/T)*tan(omegar2/2);w0=(2/T)*tan(omega0/2);B=wc2-wc1;2)将指标转换为归一化模拟低通频率指标norm_wr1=(wr12)-(w02)/(B*wr1);norm_wr2=(wr22)-(w02)/(B*wr2);norm_wc1=(wc12)-(w02)/(B*wc1);norm_wc2=(wc22)-(w02)/(B*wc2);3)设计归一化模拟低通滤波器确定归一化截频：if (abs(norm_wr1)-abs(norm_wr2)>0 norm_wr=abs(norm_wr2);else norm_wr=abs(norm_wr1);

17、endnorm_wc=1;确定归一化模拟低通滤波器阶数： N=buttord(norm_wc,norm_wr,Apmax,Armin,'s');设计归一化模拟低通滤波器： b_LP,a_LP=butter(N,norm_wc,'s');4)把归一化模拟低通滤波器转换为模拟带通滤波器b_BP,a_BP=lp2bp(b_LP,a_LP,w0,B);5)转换为数字滤波器Fs=1/T;b,a=bilinear(b_BP,a_BP,Fs);w=linspace(0,2*pi,500);h=freqz(b,a,w);plot(w,20*log10(abs(h); axis(

18、0,2*pi,-120,20); gridxlabel('frequency (rad)');ylabel('gain (dB)');sys=tf(b,a,T)例4 ：若信号由5Hz、15Hz和30Hz三个正弦频率成分构成。设计一个椭圆滤波器滤除5Hz和30Hz频率成分。解： 产生该信号Fs=100;t=(0:99)/Fs;s1=sin(2*pi*t*5); s2=sin(2*pi*t*15); s3=sin(2*pi*t*30);s=s1+s2+s3;设计该带通滤波器的通带为10Hz20Hzb,a=ellip(4,0.1,40,10 20*2/Fs);H,w=

19、freqz(b,a,512);使用FFT得到滤波前后信号的频谱S=fft(s,512);sf=filter(b,a,s);SF=fft(sf,512);作图subplot(3,2,1); plot(t,s); title('Time domain');gridylabel('original signal');subplot(3,2,5); plot(t,sf);xlabel('Time (s)');ylabe('filted signal');w1=(0:255)/256*(Fs/2);subplot(3,2,2); plot(

20、w1,abs(S(1:256); title ('Frequency domain');gridylabel('original spectrum');subplot(3,2,4); plot(w*Fs/(2*pi),abs(H);ylabel('filter');grid;subplot(3,2,6); plot(w1,abs(SF(1:256);xlabel('Frequency (Hz)');ylabel('filted spectrum');grid设计内容： 1 利用双线性变换法设计BW型及CB型数字低通

21、滤波器，绘出其幅频特性并进行比较，其指标是：（1） 通带：ApdB, Hz（2） 阻带：ArdB, =1200Hz (BW型)ArdB, =1040Hz (CB型)（3） 取样间隔T=100s2设计一个数字带通滤波器，使之满足以下指标：（1） 通带（2） =, =0.4（3） =0.8dB, 60dB2 设有一信号x(n)=1+cosn+cosn，设计各种IIR数字滤波器以达到下面的目的： 低通滤波器，滤除cosn的成分，保留的成分为1+cosn； 高通滤波器，滤除1+cosn的成分，保留的成分为cosn； 带通滤波器，滤除1+ cosn的成分，保留的成分为cosn； 带阻滤波器，滤除cosn

22、的成分，保留的成分为1+cosn； (1) 用MATLAB命令butterord定出滤波器的阶次；用butter命令设计设个滤波器；画出滤波器的幅度和相位频响；取得滤波器的系统函数H(z) （2） 试根据IIR滤波器设计原理和步骤编写出MATLAB程序。4四阶归一化低通巴特沃斯模拟滤波器系统函数为： Ha（s）用双线性变换法从Ha（s）设计四阶带通巴特沃斯数字滤波器HBP（z），并画出其幅度谱。（设T1s ，）5采用脉冲响应不变法和双线性变换法设计巴特沃斯数字低通滤波器，满足下列指标： 通带边缘频率：0.4，通带衰减：0.5dB ; 阻带边缘频率：06，阻带衰减：50dB 6用切比雪夫型设计一

23、个带通IIR数字滤波器，满足下列指标： ， ； ， ； ， 求出滤波器的系统函数H（z），并画出幅度响应图。7采用脉冲响应不变法设计一个数字低通滤波器，满足下列指标： ,通带边缘频率为1.5KHz，波动为3dB，阻带边缘频率为2KHz，衰减为40dB,通带等波纹但阻带是单调的，求出滤波器的系统函数H（z），并画出幅度响应图。8采用脉冲响应不变法设计一个数字高通滤波器，满足下列指标： ,通带边缘频率为2KHz，波动为3dB，阻带边缘频率为1.5KHz，衰减为40dB, 单调的通带和阻带，求出滤波器的系统函数H（z），并画出幅度响应图。4.2.2 窗函数的特性分析在滤波器设计中和功率谱估计中，窗函

24、数的选择对设计和分析的结果都起着很重要的作用。截短无穷长的序列会造成吉布斯现象，恰当选取窗函数，可以抑制吉布斯现象。下表给出几种常用窗函数表示式和MATLAB实现方法。 常用窗函数Matlab实现方法 函数名称MATLAB实现矩形窗(Rectangular) w=boxcar(N) 或w=ones(N,1) 汉宁窗（Hanning） w=hanning(N) 或n=0:N-1 w=1/2*(1-cos(2*pi*n/(N-1) 哈明窗(Hamming) w=hamming(N) 或 n=0:N-1 w=0.54-0.46*cos(2*pi*n/(N-1) 布莱克曼窗(Blackman) w=b

25、lackman(N) 或n=0:N-1; w=0.42-0.5*cos(2*pi*n/(N-1)+0.08*cos(4*pi*n/(N-1) Bartlett窗（三角形窗的一种） w=bartlett(N) 或 n=0:N-1; w=1-abs(2*(n-(N-1)/2)/(N-1) 凯塞窗(Kaiser) w=kaiser(N,beta) 表中前五种窗函数的形状是固定的，因而一但选择了某种窗函数，用它进行谱分析得到的频谱纹波或设计出的滤波器的阻带衰减就确定了。凯泽窗是一种可调窗，可以通过改变窗函数的形状来控制频谱纹波或阻带衰减指标，因而获得广泛的应用。 MATLAB中也提供了许多常用的窗函数

26、，调用格式为： w=boxcar(N) w= hanning(N) w = hamming(N) w = blackman(N) w =kaiser(N,beta)，N是窗函数的长度，beta是控制kaiser窗形状的参数。 例：画出布莱克曼窗；解：N=22;n=0:N-1;w=0.42-0.5*cos(2*pi*n/(N-1)+0.08*cos(4*pi*n/(N-1);% w = blackman(N)；stem(n,w);设计内容： 1分析并绘出各种窗函数的时域特性；2 使用fft函数做出各窗函数的频域特性, 并从主瓣宽度和旁瓣相对幅度两个角度进行比较分析； 3 一个序列为X(n)=co

27、s+cos，使用FFT分析其频谱： （1） 使用不同宽度的矩形窗截短该序列为M点长度, 取M分别为： a)M=20; b)M=40; c)M=160 观察不同长度的窗对谱分析结果的影响； (2 ) 使用哈明窗重做(1 ) ；对两种窗的结果进行比较和分析； 4 总结窗的不同长度和不同窗对谱分析结果的影响。 思考题： (1)在信号谱分析中，如何选择窗函数？ (2)在数字系统设计时，如何选择窗函数？ (3)如何选择不同特性的窗函数？ 4.2.3 FIR数字滤波器设计FIR数字滤波器的优点是具有线性相位。在数据通信、图象处理、语音信号处理等实际应用领域，往往都要求线性相位特性，因而使用FIR数字滤波器

28、。 FIR滤波器的传递函数表示为： H(z)=h(n) z - n其单位脉冲响应h(n)是有限长的，因此这种滤波器称为有限长单位脉冲响应(FIR)滤波器。 其设计过程的核心是求出有限的脉冲响应来逼近给定的频响。由于FIR滤波器脉冲响应在满足h(n)=±h(N-n-1)的对称的条件下具有严格的线性相位，使所设计的数字滤波器无相位失真，因此在设计中可以着重考虑幅频特性如何符合指标的要求。常用的设计方法有窗口法与频率取样法。 1窗口法设计FIR数字滤波器窗口法的设计思想是基于把给定的频响Hd(ejw)通过离散时间傅里叶逆变换(IDTFT)求得脉冲响应hd(n)，而后利用窗函数对它进行截短和

29、平滑，达到设计一个物理可实现且具有线性相位FIR数字滤波器的目的。 窗口法设计过程为： Hd(ejw)hd(n) h(n)H(ejw) Hd(ejw)是理想滤波器的频域系统函数(频响)，H(ejw)是设计出的实际滤波器的频响。加窗的作用是通过把理想滤波器的无限长脉冲响应hd(n)乘以窗函数来产生一个被截短的脉冲响应h(n)，即h(n)= hd (n) w (n)并对频响进行平滑。 例1: 分别使用矩形窗函数和哈明窗函数设计一个线性相位低通滤波器：e - j |c Hd(ej) = 0 c < | , 其中c=1rad,=12s.解: 按IDTFT求得其理想脉冲响应：hd(n) = 该理想

30、滤波器的脉冲响应hd(n)以为中心对称。根据线性相位条件h(n)= h(N-n-1)，故实际滤波器的脉冲响应长度为N=25。分别使用矩形窗函数和哈明窗函数对h(n)进行截短： 0nN-1 1) 使用矩形窗函数得h(n)= 0 其它2) 使用哈明窗函数得：0.54-0.46cos(2) 0nN-1 h(n)= 0 其它注意，滤波器阶数M=N-1=24。 %rectangular window加矩形窗 M=24;n=0:M;omegac=1; %M: order of the filter hd=omegac/pi.*sinc(omegac/pi*(n-M/2); stem(0:M,hd); h1

31、=hd; H1=fft(h1,512); H1_db=20*log10(abs(H1); %hamming window,哈明窗 w=0.54-0.46.*cos(2*pi*n./M); stem(0:M,w); 加哈明窗 h2=hd.*w; H2=fft(h2,512); H2_db=20*log10(abs(H2); 做FIR滤波器的幅频特性曲线 omega=0:pi/255:pi; c=plot(omega,H1_db(1:256),omega,H2_db(1:256),r); 上图绘出了使用矩形窗、哈明窗设计的实际滤波器的幅频特性。 使用信号处理工具箱中的fir1函数可以很容易地实现F

32、IR滤波器窗口法设计。它可设计低通、高通、带通、带阻滤波器及多频带滤波器。格式如下： b = fir1(n,Wc) b = fir1(n,Wc,ftype) b = fir1(n,Wc,window) b = fir1(n,Wc,ftype,window) b = fir1(n,Wc,ftype,window,noscale) n为滤波器的阶数，Wc为3dB截频：0 < Wc < 1, 1对应取样频率的一半。对带通和带阻滤波器，Wc=w1,w2，为通带上下截频。若Wc=w1,w2，w3wn,表示滤波器具有n个频带：0<w<w1, w2<w<w3,wn<

33、w<1。 ftype指定滤波器类型： ftype=high, 指定一个截频为Wc的高通滤波器； ftype=stop 指定一个带阻滤波器，其阻带截止频率为Wc=w1,w2； ftype=DC-0 在多频带滤波器中，使第一个频带0<w<w1为阻带； ftype=DC-1 在多频带滤波器中，使第一个频带0<w<w1为通带； window为指定的窗向量，若不指定，默认为哈明窗。 使用信号处理工具箱中的fir1函数对本例题进行设计： %rectangular window b1=fir1(24,1/pi,boxcar(25); H1,w=freqz(b1,1,512);

34、H1_db=20*log10(abs(H1); %hamming window b2=fir1(24,1/pi,hamming(25); H2,w=freqz(b2,1,512); H2_db=20*log10(abs(H2); c=plot(w,H1_db,w,H2_db,r); 2. 频率取样法FIR滤波器 频率取样法的设计思想是基于对给出的理想频响Hd(ejw) 进行取样，通过离散傅里叶反变换(IDFT)从频谱样点直接求得有限脉冲响应。 频率取样法的过程如下： Hd(ejw) H(k) h(n) H(ejw) 其中Hd(ejw)是理想滤波器的数字域系统函数(频响)，H(ejw)是设计出来

35、的实际滤波器系统函数。频域取样法的关键是正确确定频域取样点，在0-2内的样点有如下约束条件： H(k)=H(N-k) (k)= - (N-k)例2： 采用频率取样法设计FIR滤波器, 使之满足 e - j , |cHd(ej) = 0 , c<| 其中，c = 1rad, = 12s.解：N为奇数且h (n)为偶对称，故采用N为奇数的I型取样。设计程序如下：N=25; tao=(N-1)/2; omegac=1;m=fix(omegac/(2*pi/N)+1);omega=0:N-1*2*pi/N;ABS_H =ones(1,m), zeros(1,N-2*m+1),ones(1,m-1

36、);H=ABS_H.*exp(-j*tao.*omega);h=ifft(H);figure(1)stem(omega,abs(H); title('FIR filter')ylabel('mag of Frequency response'); xlabel('Frequency (rad)')hold on;i=0:511*2*pi/512; m1=fix(omegac/(2*pi/512)+1);ABS_HD=ones(1,m1),zeros(1,512-2*m1+1),ones(1,m1-1);c=plot(i,ABS_HD,'r

37、'); legend(c,'ideal mag.',0);hold off;figure(2)stem(0:N-1,real(h);title('impulse response'); xlabel('Tume n'); ylabel('hn');Hh=fft(h, 512)Hh_db=20*log10(abs(Hh);figure(3)omega=0:255*pi/256;plot(omega,Hh_db(1:256); gridset(gca, 'ytick', (-120:20:-20),-3,0);xlabel('Frequency (rad)'); ylabel('Norm