数字滤波器设计

24第四章数字滤波其次节数字滤波的设计通信与电子信息当中，在对信号作分析与处理时，常会用到有用信号叠加无用噪声的问题。这些噪声信号有的是与信号同时产生的，有的是在传输过程中混入的，在接收的信号中，必需消退或减弱噪声干扰，这是信号处理中格外重要的问题。根据有用信号与噪声的不同特性，消退或减弱噪声，提取有用信号的过程就称为滤波。滤波器的种类很多，实现方法也多种多样，本章利用Matlab来进展数字滤波器的设计。数字滤波器是一离散时间系统，它对输入序列x(n)进展加工处理后，输出序列y(n)y(n)x(n)的频谱相比发生某种变化。DSP理论得知，无限长冲激响应〔IIR〕需要递归模型来实现，有限长冲激响应〔FIR〕滤波器可以承受递归的方式也可承受非递归的方式实现。本章把FIRIIR滤波器分别用Matlab进展分析与设计。数字滤波器的构造参看《数字信号处理》一书。数字滤波器的设计一般经过三个步骤：给出所需滤波器的技术指标。设计一个H〔，使其靠近所需要的技术指标。实现所设计的H〔。IIRp IIR数字滤波器的任务就是寻求一个因果、物理可实现的系统函数H(z)，使它的频响H(ejω)满足所期望得到的低通频域指标，即通带衰减A、阻带衰减Ap 通带截频ω、阻带截频ω。而其它形式的滤波器由低通的变化得到。p r承受间接法设计IIR数字滤波器就是按给定的指标，先设计一个模拟滤波器，进而通过模拟域与数字域的变换，求得物理可实现的数字滤波器。从模拟滤波器变换到数字滤波器常用的有：脉冲响应不变法和双线性变换法。IIR滤波器的设计过程如下数字频域指标模拟频域指标设计模拟滤波器H(S) 设计数字滤波器H(z)1.模拟滤波器简介模拟滤波器的设计方法已经进展得格外成熟，常用的高性能模拟低通滤波器有巴特沃斯型、切比雪夫型和椭圆型，而高通、带通、带阻滤波器则可以通过对低通进展频率变换来求得。必需指出，这三种滤波器都是非线性的相频特性。具体设计方法参见《数字信号处理》一书。巴特沃斯低通滤波器的特点是：通、阻带均为单调下降。这种单调下降的特性使得系统的误差分布不均匀。在设计中，假设在通带满足指标，阻带指标就过于好了。或者说，在阻带满足指标，通带指标也过于好了。这样不利于以最小的阶数来满足设计指标。切比雪夫模拟低通滤波器分成Ι型和ΙΙ型：Ι型是通带等波浪，阻带单调下降；ΙΙ型滤波器是通带单调下降，阻带等波浪。切比雪夫的等波浪特性使得可以用较小的阶数设计出符合指标的滤波器。椭圆滤波器特点是通、阻带均为等波浪。椭圆滤波器的误差均匀分布在通、阻带，比较上面几种滤波器，在同样衰耗指标下，椭圆滤波器所要求的阶数最小。在同样阶数状况下，椭圆滤波器的通带到阻带的变化最陡峭、性能最好。2．模拟滤波器映射到数字滤波器从模拟滤波器变换到数字滤波器常用的有：冲激响应不变法和双线性变换法。〔1〕冲激响应不变法aa冲激响应不变法的根本思路为：设模拟滤波器的单位冲激响应为h〔t〕，数h(n)h〔t〕的均匀取样，即aaH(S)逆换h(t)样 h(n)=h(nT) 换H(z)冲激响应不变法由于保持了模拟滤波器的冲激响应在取样时刻的取样值，所以具有较好的时间特性。但由于从模拟域〔s域〕变换到数字域〔z域〕的映射关系不是一一对应的关系，因而频谱混叠难以避开,所以在实际应用中有它的局限性。3．Matlab滤波器设计函数简介MATLAB信号处理工具箱中供给的函数可以很简洁地实现IIR滤波器设计。常用的IIR数字滤波器设计函数有：IR滤波器阶数〔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’,’s’,分别对应低通、高通、带通、带阻、模拟滤波器，默认状况下，为低通或带通。对于带通和带阻滤波器Wp=[Wp1,Wp2]，Wp1和Wp2为通带的上、下截频，Wr=[Wr1,Wr2]，Wr1Wr2分别为阻带的上、Options=’s’,则Wp和Wr以〔归一或非归一〕1。输出参数N为滤波器的阶数，Wc3dB截频。假设Ap=3dB，则Wc=Wp。其他滤波器的定阶函数类似buttord的格式：切比雪夫I型滤波器[N,Wc]=cheb1ord(Wp,Wr,Ap,Ar)切比雪夫II型滤波器[N,Wc]=cheb2ord(Wp,Wr,Ap,Ar)椭圆滤波器[N,Wc]=ellipord(Wp,Wr,Ap,Ar)完整的滤波器设计函数，依据滤波器阶数和通带截频设计出实际数字滤波器H(z)〔模拟域到数字域的变换只承受双线性变换法〕：butter-巴特沃斯〔Butterworth〕滤波器设计cheby1-切比雪夫〔Chebyshev〕I型滤波器设计cheby2--切比雪夫〔Chebyshev〕II型滤波器设计ellip-椭圆〔Elliptic〕滤波器设计maxflat-通用的巴特沃斯〔Butterworth〕低通滤波器设计yulewalk-Yule-Walker滤波器设计〔直接数字滤波器设计法〕butter的格式如下：[b,a]=butter(N,Wc,options)[z,p,k]=butter(N,Wc,options)[A,B,C,D]=butter(N,Wc,options)N为滤波器的阶数，Wc3dB截频：0<Wc<1。Options=’low’,’high’,’bandpass’,’stop’,’s’,分别对应低通、高通、带通、带阻、模拟滤波器，默认状况下，为低通或带通。对于带通和带阻滤波器，Wc=[W1,W2]，W1和W2分别为通带的上、下截频。输出变量的三种格式对应着滤波器系统函数H(z)的分子分母多项式形式、零极点-增益形式和状态空间形式。假设设计模拟滤波器，令Options=’s’,则Wc以rad/s〔归一化或非归一化〕给出，因此可大于1。例1．0~0.2π3dB0.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.5π02〕 上截频 =0.55π,下截频 ＝0.45πc2 c13〕 在阻带r ＝0.6π处，最小衰耗A＝10dB设取样周期＝5π1s要求｜｜>衰耗曲线单调增加。0解：〔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=omegar2/T;B=wc2-wc1;将指标转换为归一化模拟低通频率指标norm_wr1=(((wr1^2)-(w0^2))/(B*wr1));norm_wr2=(((wr2^2)-(w0^2))/(B*wr2));norm_wc1=(((wc1^2)-(w0^2))/(B*wc1));norm_wc2=(((wc2^2)-(w0^2))/(B*wc2));设计归一化模拟低通滤波器确定归一化截频if(abs(norm_wr1)-abs(norm_wr2))>0norm_wr=abs(norm_wr2);elsenorm_wr=abs(norm_wr1);endnorm_wc=1;确定归一化模拟低通滤波器阶数N=buttord(norm_wc,norm_wr,Apmax,Armin,”s”);设计归一化模拟低通滤波器[b_LP,a_LP]=butter(N,norm_wc,”s”);把归一化模拟低通滤波器转换为模拟带通滤波器[b_BP,a_BP]=lp2bp(b_LP,a_LP,w0,B);利用脉冲响应不变法把模拟滤波器转换为数字滤波器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]);gridxlabel(”frequency(rad)”);ylabel(”gain(dB)”);sys=tf(b,a,T)(2)双线性变换法双线性变换法抑制了脉冲响应不变法存在频谱混叠的缺点,但在设计过程中确定要留意，由于宽频带压缩到较窄的频率范围，导致了频率畸变。其根本思想是：首先按给定的指标设计一个模拟滤波器，通过适当的数学变换方法将s域映射到z域，把无限宽的频带变换成有限宽频带。由于在实现数字化以前已经对频带进展了压缩，所以在数字化以后的频响可以做到无频谱混叠。该法可以用来设计各种不同频带要求的数字滤波器，所以获得广泛应用。例3：按[2]的指标使用双线性变换法设计一个数字带通滤波器。解：〔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=(2/T)*tan(omegar2/2);w0=(2/T)*tan(omega0/2);B=wc2-wc1;将指标转换为归一化模拟低通频率指标norm_wr1=(((wr1^2)-(w0^2))/(B*wr1));norm_wr2=(((wr2^2)-(w0^2))/(B*wr2));norm_wc1=(((wc1^2)-(w0^2))/(B*wc1));norm_wc2=(((wc2^2)-(w0^2))/(B*wc2));设计归一化模拟低通滤波器确定归一化截频：if(abs(norm_wr1)-abs(norm_wr2))>0norm_wr=abs(norm_wr2);elsenorm_wr=abs(norm_wr1);endnorm_wc=1;确定归一化模拟低通滤波器阶数：N=buttord(norm_wc,norm_wr,Apmax,Armin,”s”);设计归一化模拟低通滤波器：[b_LP,a_LP]=butter(N,norm_wc,”s”);把归一化模拟低通滤波器转换为模拟带通滤波器[b_BP,a_BP]=lp2bp(b_LP,a_LP,w0,B);转换为数字滤波器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([0,2*pi,-120,20]);gridxlabel(”frequency(rad)”);ylabel(”gain(dB)”);sys=tf(b,a,T)45Hz、15Hz30Hz三个正弦频率成分构成。设计一个椭圆滤波器5Hz30Hz频率成分。解：产生该信号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;设计该带通滤波器的通带为10Hz－20Hz[b,a]=ellip(4,0.1,40,[1020]*2/Fs);[H,w]=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(”Timedomain”);gridylabel(”originalsignal”);subplot(3,2,5); plot(t,sf);xlabel(”Time (s)”);ylabe(”filtedsignal”);w1=(0:255)/256*(Fs/2);subplot(3,2,2); plot(w1,abs(S(1:256)));title(”Frequencydomain”);gridylabel(”originalspectrum”);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(”filtedspectrum”);grid：利用双线性变换法设计BW型及CB型数字低通滤波器，绘出其幅频特性并进展比较，其指标是：〔1〕通带：Ap0.7dB, 01000Hz〔2〕阻带：Ar10dB, r=1200Hz (BW型)Ar10dB, =1040Hz (CB型)r取样间隔T=100s设计一个数字带通滤波器，使之满足以下指标：〔1〕通带0.250.35〔2〕 , =0.4l u〔3〕 =0.8dB, 60dBp r2x(n)=1+cosn+cos2

n，设计各种IIR数字滤波器以到达下面的4 3目的：2 低通滤波器，滤除cos3n1+cos4n；1+cos

24n的成分，保存的成分为cos3n；2 1+cos

3n的成分，保存的成分为cos4n；带阻滤波器，滤除cos

24n1+cos3n；MATLAB命令butterord定出滤波器的阶次；用butter命令设计设个滤波器；画出滤波器的幅度和相位频响；取得滤波器的系统函数H(z)〔2〕试依据IIR滤波器设计原理和步骤编写出MATLAB程序。四阶归一化低通巴特沃斯模拟滤波器系统函数为：1Ha〔s〕＝

s42.6131s33.4142s22.6131s1用双线性变换法从Ha〔s〕设计四阶带通巴特沃斯数字滤波器HBP〔z〕，并画出其幅度谱。〔设T＝1s1c

0.35，2c

0.65〕承受脉冲响应不变法和双线性变换法设计巴特沃斯数字低通滤波器，满足以下指标：通带边缘频率：0.4，通带衰减：0.5dB ;阻带边缘频率：06，阻带衰减：50dB用切比雪夫型设计一个带通IIR数字滤波器，满足以下指标：|0.25 ；0|H(ej)|0.01 ；0.95|H(ej)|1.05，0.75||求出滤波器的系统函数H〔z〕，并画出幅度响应图。7．承受脉冲响应不变法设计一个数字低通滤波器，满足以下指标：f=8KHz,1.5KHz3dB2KHz，衰减s为40dB,通带等波浪但阻带是单调的，求出滤波器的系统函数H〔z〕，并画出幅度响应图。8．承受脉冲响应不变法设计一个数字高通滤波器，满足以下指标：f=1KHz,2KHz3dB1.5KHz，衰减s40dB,单调的通带和阻带，求出滤波器的系统函数H〔z〕，并画出幅度响应图。窗函数的特性分析在滤波器设计中和功率谱估量中，窗函数的选择对设计和分析的结果都起着很重要的作用。截短无穷长的序列会造成吉布斯现象，恰中选取窗函数，可以抑制吉布斯现象。下表给出几种常用窗函数表示式和MATLAB实现方法。常用窗函数Matlab实现方法函数名称函数名称MATLAB实现矩形窗(Rectangular)w=boxcar(N)或w=ones(N,1)汉宁窗〔Hanning〕w=hanning(N)或n=0:N-1w=1/2*(1-cos(2*pi*n/(N-1)))哈明窗(Hamming)w=hamming(N)或n=0:N-1w=0.54-0.46*cos(2*pi*n/(N-1))布莱克曼窗(Blackman)w=blackman(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))w=1-abs(2*(n-(N-1)/2)/(N-1))凯塞窗(Kaiser)w=kaiser(N,beta)表中前五种窗函数的外形是固定的，因而一但选择了某种窗函数，用它进展谱分析得到的频谱纹波或设计出的滤波器的阻带衰减就确定了。凯泽窗是一种可调窗，可以通过转变窗函数的外形来把握频谱纹波或阻带衰减指标，因而获得广泛的应用。MATLAB中也供给了很多常用的窗函数，调用格式为：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);设计内容：分析并绘出各种窗函数的时域特性；fft函数做出各窗函数的频域特性,并从主瓣宽度和旁瓣相对幅度两个角度进展比较分析；11一个序列为X(n)=cos

n，使用FFT分析其频谱：20 20〔1〕使用不同宽度的矩形窗截短该序列为M点长度,取M分别为：a)M=20; b)M=40; c)M=160观看不同长度的窗对谱分析结果的影响；(2)使用哈明窗重做(1)；对两种窗的结果进展比较和分析；总结窗的不同长度和不同窗对谱分析结果的影响。思考题：在信号谱分析中，如何选择窗函数？在数字系统设计时，如何选择窗函数？如何选择不同特性的窗函数？FIRFIR数字滤波器的优点是具有线性相位。在数据通信、图象处理、语音信号处理等实际应用领域，往往都要求线性相位特性，因而使用FIR数字滤波器。FIR滤波器的传递函数表示为：H(z)=N1 h(n)z-nn0其单位脉冲响应h(n)是有限长的，因此这种滤波器称为有限长单位脉冲响应(FIR)滤波器。其设计过程的核心是求出有限的脉冲响应来靠近给定的频响。由于FIR滤波器脉冲响应在满足h(n)=±h(N-n-1)的对称的条件下具有严格的线性相位，使所设计的数字滤波器无相位失真，因此在设计中可以着重考虑幅频特性如何符合指标的要求。常用的设计方法有窗口法与频率取样法。窗口法设计FIR数字滤波器窗口法的设计思想是基于把给定的频响Hd(ejw)通过离散时间傅里叶逆变换(IDTFT)求得脉冲响应hd(n)，而后利用窗函数对它进展截短和平滑，到达设计一个物理可实现且具有线性相位FIR数字滤波器的目的。窗口法设计过程为：d H(j)DFTh(n) ()h(n)FTH(jd Hd(ejw)是抱负滤波器的频域系统函数(频响)，H(ejw)是设计出的实际滤波器的频响。加窗的作用是通过把抱负滤波器的无限长脉冲响应hd(n)乘以窗函数来产生一个被h(n)h(n)=hd(n)w(n)并对频响进展平滑。1:分别使用矩形窗函数和哈明窗函数设计一个线性相位低通滤波器：c-jω |ω|ωcHd(ejω)=

0 ωc<|ω|

, 其中ωc=1rad,=12s.IDTFT求得其抱负脉冲响应：h

csinc( c(n))d d该抱负滤波器的脉冲响应h(n)以为中心对称。依据线性相位条件h(n)=h(N-n-1)，故实际滤波器的脉冲响应长度为N=25。分别使用矩形窗函数和哈明窗函数对h(n)进展截短：dcsinc(c

c(n)) 0nN-1h(n)=0 其它使用哈明窗函数得：ccsinc(

c(n))[0.54-0.46cos(2

)] 0nN-1N1h(n)=0 其它留意，滤波器阶数M=N-1=24。%rectangularwindow加矩形窗M=24;n=0:M;omegac=1; %M:orderofthefilterhd=omegac/pi.*sinc(omegac/pi*(n-M/2));stem([0:M],hd);h1=hd;H1=fft(h1,512);H1_db=20*log10(abs(H1));%hammingwindow,哈明窗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函数可以很简洁地实现FIR滤波器窗口法设计。它可设计低通、高通、带通、带阻滤波器及多频带滤波器。格式如下：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为滤波器的阶数，Wc3dB截频：0<Wc<1,1对应取样频率的一半。对带通和带阻滤波器，Wc=[w1,w2]，为通带上下截频。假设Wc=[w1,w2，w3…wn],表示滤波器具有n个频带：0<w<w1,w2<w<w3,…wn<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函数对本例题进展设计：%rectangularwindowb1=fir1(24,1/pi,boxcar(25));[H1,w]=freqz(b1,1,512);H1_db=20*log10(abs(H1));%hammingwindowb2=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—‘);频率取样法FIR滤波器频率取样法的设计思想是基于对给出的抱负频响Hd(ejw)进展取样，通过离散傅里叶反变换(IDFT)从频谱样点直接求得有限脉冲响应。频率取样法的过程如下：频域取样(频域取样(2K)dH(j) N H(k) FTh(n) FTH(ej)dHd(ejw)是抱负滤波器的数字域系统函数(频响)，H(ejw)是设计出来的实际滤波器系统函数。频域取样法的关键是正确确定频域取样点，在ω∈[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)];H=ABS_H.*exp(-j*tao.*omega);h=ifft(H);figure(1)stem(omega,abs(H));title(”FIRfilter”)ylabel(”magofFrequencyresponse”);xlabel(”Frequency(rad)”)holdon;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”);legend(c,”idealmag.”,0);holdoff;figure(2)stem([0:N-1],real(h));title(”impulseresponse”);xlabel(”Tumen”);ylabel(”h[n]”);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(”Normalizedmagnitude(dB)”);title(”Frequencyresponse”);例３：承受I型频率取样法设计FIR滤波器,满