精品资料欢迎下载实验四窗函数法设计FIR 数字滤波器一、实验目的1、掌握窗函数法设计FIR 数字滤波器的原理及具体方法。 2、掌握频率取样法设计FIR 数字滤波器的原理和基本方法。 3、学习利用窗函数法和频率取样法设计低通、带通、高通、带阻数字滤波器。二、实验环境计算机、 MATLAB 软件三、实验基础理论窗函数设计FIR 滤波器1. 基本原理窗函数设计法的基本思想为，首先选择一个适当的理想的滤波器()jdHe，然后用窗函数截取它的单位脉冲响应(n)dh，得到线性相位和因果的FIR 滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器，使设计的滤波器的单位脉冲响应逼近理想滤波器的单位脉冲响应。2. 设计步骤（1）给定理想滤波器的频率响应()jdHe，在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽为()cc的低通滤波器由下式给定：|,0)(,| ,)(cjdcjajdeHeeH其中为采样延迟，其作用是为了得到一个因果系统。（2）确定这个滤波器的单位脉冲响应)()(sin()(anannhcd为了得到一个(n)h长度为 N的因果的线性相位FIR 滤波器，我们令21Na（3）用窗函数截取(n)dh得到所设计FIR 数字滤波器：)()()(nRnhnhNd3. 窗函数的选择常用的窗函数有矩形（Rectangular ）窗，汉宁（Hanning ）窗，海明（ Hamming ）窗、布莱克曼（ Blackman）窗、凯瑟（Kaiser ）窗等表 4-1 MATLAB中产生窗函数的命令MATLAB 函数窗函数MATLAB 函数窗函数Boxcar 矩形窗函数Blackman 布莱克曼窗Hanning 汉宁窗函数Kaiser 凯瑟窗函数Hamming 海明窗精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 16 页精品资料欢迎下载表 4-2 常用窗函数的特性窗函数窗函数频率特性加窗后滤波器指标旁瓣峰值dB 主瓣宽度过渡带宽最小阻带衰减dB 矩形窗-13 4/N 1.8 /N -21 汉宁窗-31 8/N 6.2 /N -44 海明窗-41 8/N 6.6 /N -53 布莱克曼窗-57 12/N 11/N -74 凯瑟窗是一种广泛在实际中广泛应用的窗函数，它由下式给定：20021(1)1( )nINnI其中0Ix是修正的零阶贝塞尔函数，参数控制最小阻带衰减，这种窗函数对于相同的 N可以提供不同的过渡带宽。由于贝塞尔函数比较复杂，这种窗函数的设计方程很难推导，然而幸运的是，有一些经验设计方程可以直接使用。已知给定的指标,RpstpsA和，滤波器长度N和凯瑟窗参数可以按如下凯瑟窗方程给出过渡带带宽：stp7.9512.285sAN0.40.1102(8.7),500.5842(21)0.07886(21),2150sssssAAAAA频率取样设计FIR 滤波器1. 基本原理频率取样法从频域出发，把理想的滤波器()jdHe等间隔采样得到( )dHk, 将( )dHk作为实际设计滤波器的( )H k：2( )( )() |0,1,1jdkNH kHkH ekN得到( )H k以后可以由( )H k来确定唯一确定滤波器的单位脉冲响应( )h n,()jH e可以由( )H k求得：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 16 页精品资料欢迎下载10( )( )2()( ) ()Njkh nIDFT H kH eH kkN其中( )x为内插函数：12sin(/ 2)sin(/ 2)NjNeN（有( )H k求得的频率响应()jH e将逼近()jdHe。如果我们设计的是线性相位FIR 滤波器， 则( )H k的幅度和相位满足线性相位滤波器的约束条件。我们将( )H k表示为如下形式( )( )( )=( )( )jkjkrH kH k eHk e当( )h n为实数，则*( )()H kHNk由此得到( )()rrHkHNk即( )/ 2rHkkN以为中心偶对称。在利用线性相位条件可知，对于1 型和 2 型线性相位滤波器：1 21()0,22( )1 21()()1,122NkNkNkNNNkkNN对于 3 型和 4 型线性相位滤波器1 21()0,222( )1 21()()1,1222NkNkNkNNNkkNN2. 设计步骤（1）由给定的理想滤波器给出( )rHk和( )k。（2）由()( )( )=( )( )jkjkrH kH k eHk e求得( )H k精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 16 页精品资料欢迎下载（3）根据( )H k求得( )h n或()jH e四、实验内容1、设计一个数字低通FIR 滤波器，其技术指标如下：0.2 ,0.25ppRdB0.3 ,50stsAdB分别采用矩形窗、 汉宁窗、 海明窗、布莱克曼窗、 凯瑟窗设计该滤波器。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。（1）矩形窗程序代码：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp; N=ceil(1.8*pi/tr_width) n=0:N-1; wc=(wst+wp)/2; alpha=(N-1)/2; hd=(wc/pi)*sinc(wc/pi)*(n-alpha); w_boxcar=boxcar(N); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n );ylabel(hd(n); Hr,w1=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(omega/pi);ylabel(H(omega); subplot(223); stem(n,h,filled); axis tight;xlabel(n );ylabel(h(n); H,w=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H); axis tight;xlabel(omega/pi);ylabel(dB ); grid on; MATLAB 图形：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 16 页精品资料欢迎下载（2）汉宁窗程序代码：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp; N=ceil(6.2*pi/tr_width) n=0:N-1; wc=(wst+wp)/2; alpha=(N-1)/2; hd=(wc/pi)*sinc(wc/pi)*(n-alpha); w_boxcar=hanning(N); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n );ylabel(hd(n); Hr,w1=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(omega/pi);ylabel(H(omega); subplot(223); stem(n,h,filled); axis tight;xlabel(n );ylabel(h(n); H,w=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H); 05101500.10.2nhd(n)00.51-0.500.511.5/H()05101500.10.2nh(n)00.20.40.60.8-80-60-40-200/dB精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 16 页精品资料欢迎下载axis tight;xlabel(omega/pi);ylabel(dB ); grid on; MATLAB 图形：（3）海明窗程序代码：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp; N=ceil(6.6*pi/tr_width) n=0:N-1; wc=(wst+wp)/2;alpha=(N-1)/2; hd=(wc/pi)*sinc(wc/pi)*(n-alpha); w_boxcar=hamming(N); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n );ylabel(hd(n); Hr,w1=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(omega/pi);ylabel(H(omega); subplot(223); stem(n,h,filled); axis tight; xlabel(n ); 020406000.10.2nhd(n)00.51-0.500.511.5/H()020406000.050.10.150.2nh(n)00.20.40.60.8-120-100-80-60-40-200/dB精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 16 页精品资料欢迎下载ylabel(h(n); H,w=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H); axis tight;xlabel(omega/pi);ylabel(dB ); grid on; MATLAB 图形：（4）布莱克曼窗程序代码：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp; N=ceil(11*pi/tr_width) n=0:N-1; wc=(wst+wp)/2;alpha=(N-1)/2; hd=(wc/pi)*sinc(wc/pi)*(n-alpha); w_boxcar=blackman(N); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n );ylabel(hd(n); Hr,w1=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(omega/pi);ylabel(H(omega); subplot(223); stem(n,h,filled); axis tight;xlabel(n );ylabel(h(n); 020406000.10.2nhd(n)00.51-0.500.511.5/H()020406000.10.2nh(n)00.20.40.60.8-100-500/dB精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 16 页精品资料欢迎下载H,w=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H); axis tight;xlabel(omega/pi);ylabel(dB ); grid on; MATLAB 图形为：（5）凯瑟窗程序代码：wp=0.2*pi;wst=0.3*pi;tr_width=wst-wp;As=50; N=ceil(As-7.95)/(2.285*tr_width)+1; beta=0.1102*(As-8.7); n=0:N-1; wc=(wst+wp)/2;alpha=(N-1)/2; hd=(wc/pi)*sinc(wc/pi)*(n-alpha); w_boxcar=kaiser(N,beta); h=hd.*w_boxcar; subplot(221); stem(n,hd,filled); axis tight;xlabel(n );ylabel(hd(n); Hr,w1=zerophase(h); subplot(222); plot(w1/pi,Hr); axis;xlabel(omega/pi);ylabel(H(omega); subplot(223); stem(n,h,filled); 05010000.10.2nhd(n)00.51-0.500.511.5/H()05010000.10.2nh(n)00.20.40.60.8-150-100-500/dB精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 16 页精品资料欢迎下载axis tight;xlabel(n );ylabel(h(n); H,w=freqz(h,1); subplot(224); plot(w/pi,20*log10(abs(H)/max(H); axis tight;xlabel(omega/pi);ylabel(dB ); grid on; MATLAB 图形：2、设计一个数字带通FIR 滤波器，其技术指标如下：下阻带边缘：10.2 ,60stsAdB下通带边缘：10.35 ,1ppRdB上通带边缘：20.65 ,1ppRdB上阻带边缘：20.8 ,60stsAdB程序代码：wp1=0.2*pi;Rp1=1; wst1=0.35*pi;A1=60; width1=wst1-wp1; N1=ceil(11*pi/width1)+1; n1=0:(N1-1); wc1=(wp1+wst1)/2; alpha=(N1-1)/2; wp2=0.65*pi;Rp2=1;wst2=0.8*pi;A2=60; 0204000.10.2nhd(n)00.51-0.500.511.5/H()0204000.050.10.150.2nh(n)00.20.40.60.8-100-500/dB精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 16 页精品资料欢迎下载width2=wst2-wp2; N2=ceil(11*pi/width2)+1; n2=0:(N2-1); wc2=(wp2+wst2)/2; alpha=(N2-1)/2; hd=(wc2/pi)*sinc(wc2/pi)*(n2-alpha)-(wc1/pi)*sinc(wc1/pi)*(n1-alpha); w_w=blackman(N1); h=hd.*w_w; subplot(221); stem(n1,h,filled); subplot(222); H,w=freqz(h,1); plot(w/pi,20*log10(abs(H)/max(abs(H); subplot(223); Hr,w1=zerophase(h); plot(w1/pi,Hr); subplot(224); stem(n1,hd,filled); Hr,wl=zerophase(h); grid on; MATLAB 图形： 3. 采用频率取样法设计FIR 数字低通滤波器，满足以下指标020406080-0.4-0.200.20.40.600.51-150-100-50000.51-0.500.511.5020406080-0.4-0.200.20.40.6精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 16 页精品资料欢迎下载0.2 ,0.250.3 ,50ppstsRdBAdB（1）取 N=20，过渡带没有样本。（2）取 N=40，过渡带有一个样本，T=0.39 。（3）取 N=60，过渡带有两个样本，T1=0.5925 ，T2=0.1099 。（4）分别讨论采用上述方法设计的数字滤波器是否都能满足给定的指标要求。（1）程序代码：N=20; alpha=(N-1)/2; L=0:N-1;wL=(2*pi/N)*L; Hrs=1,1,1,zeros(1,15),1,1; Hdr=1,1,0,0; wdL=0,0.25,0.25,1; k1=0:floor(N-1)/2); k2=floor(N-1)/2)+1:N-1; angH=-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2); H=Hrs.*exp(j*angH); h=ifft(H,N); w=0:500*pi/500; Hr,wr=zerophase(h); subplot(221); plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1); axis(0,1,-0.1,1.1); xlabel(omega (pi); ylabel(Hr(k); subplot(222); stem(L,h,filled); axis(0,N-1,-0.1,0.3); xlabel(n ); ylabel(h(n); subplot(223); plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1); axis(0,1,-0.2,1.2); grid on; xlabel(omega (pi); ylabel(Hr(omega); subplot(224); plot(wr/pi,20*log10(abs(Hr)/max(abs(Hr); axis(0,1,-50,5); grid on;xlabel(omega (pi);ylabel(dB ) MATLAB 图形如下：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 16 页精品资料欢迎下载（2）程序代码：N=40; alpha=(N-1)/2; L=0:N-1; wL=(2*pi/N)*L; Hrs=1,1,1,1,1,0.39,zeros(1,29),0.39,1,1,1,1; Hdr=1,1,0.39,0,0; wdL=0,0.2,0.25,0.3,1; k1=0:floor(N-1)/2); k2=floor(N-1)/2)+1:N-1; angH=-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2); H=Hrs.*exp(j*angH); h=ifft(H,N); w=0:500*pi/500; Hr,wr=zerophase(h); subplot(221); plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1); axis(0,1,-0.1,1.1); xlabel(omega (pi); ylabel(Hr(k); subplot(222); stem(L,h,filled); axis(0,N-1,-0.1,0.3); xlabel(n ); 00.5100.51 ( )Hr(k)051015-0.100.10.20.3nh(n)00.51-40-200 ( )dB00.5100.51 ( )Hr()精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 16 页精品资料欢迎下载ylabel(h(n); subplot(223); plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1); axis(0,1,-0.2,1.2); grid on; xlabel(omega (pi); ylabel(Hr(omega); subplot(224); plot(wr/pi,20*log10(abs(Hr)/max(abs(Hr); axis(0,1,-50,5); grid on; xlabel(omega (pi); ylabel(dB ) grid on; MATLAB 图形如下：3、程序代码：N=60; alpha=(N-1)/2; L=0:N-1; wL=(2*pi/N)*L; Hrs=1,1,1,1,1,1,1,0.5925,0.1099,zeros(1,43),0.1099,0.5925,1,1,1,1,1,1; Hdr=1,1,0.5925,0.1099,0,0; wdL=0,0.2,7/30,8/30,0.3,1; k1=0:floor(N-1)/2); 00.5100.51 ( )Hr(k)0102030-0.100.10.20.3nh(n)00.5100.51 ( )Hr()00.51-40-200 ( )dB精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 16 页精品资料欢迎下载k2=floor(N-1)/2)+1:N-1; angH=-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2); H=Hrs.*exp(j*angH); h=ifft(H,N);w=0:500*pi/500; Hr,wr=zerophase(h); subplot(221); plot(wdL,Hdr,wL(1:N/2+1)/pi,Hrs(N/2+1); axis(0,1,-0.1,1.1); xlabel(omega (pi); ylabel(Hr(k); subplot(222); stem(L,h,filled); axis(0,N-1,-0.1,0.3); xlabel(n ); ylabel(h(n); subplot(223); plot(wr/pi,Hr,wL(1:N/2+1)/pi,Hrs(1:N/2+1); axis(0,1,-0.2,1.2); grid on; xlabel(omega (pi); ylabel(Hr(omega); subplot(224); plot(wr/pi,20*log10(abs(Hr)/max(abs(Hr); axis(0,1,-50,5);grid on;xlabel(omega (pi);ylabel(dB ) grid on; MATLAB 图形 : 00.5100.51 ()Hr(k)02040-0.100.10.20.3nh(n)00.5100.51 ()Hr()00.51-40-200 ( )dB精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 16 页精品资料欢迎下载各阶数的通阻带指标如下表：阶数通带波动 0.25dB 阻带衰减50dB 是否满足指标N=20 不满足不满足否N=40 不满足满足否N=60 满足满足是4. 采用频率取样技术设计下面的高通滤波器0.6,500.8,1stsppAdBRdB对于高通滤波器，N必须为奇数（或1 型滤波器）。选择 N=33，过渡带有两个样本，过渡带的最优值为T1=0.1095 ，T2=0.598. 程序代码：N=33;alpha=(N-1)/2;L=0:N-1;wL=(2*pi/N)*L; Hrs=0,0,0,0,0,0,0,0,0,0,0,0.1095,0.598,1,1,1,1,1,1,1,0.598,0.1095,0,0,0,0,0,0,0,0,0,0,0; Hdr=0,0,0.1095,0.598,1,1; wdL=0,0.6,2/3,24/33,0.8,1; k1=0:floor(N-1)/2); k2=floor(N-1)/2)+1:N-1; angH=-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2); H=Hrs.*exp(j*angH); h=ifft(H,N); w=0:500*pi/500; Hr=h*cos(alpha-L)*w); subplot(221); plot(wdL,Hdr,wL(1:17)/pi,Hrs(1:17); axis(0,1,-0.1,1.1); xlabel(omega (pi); ylabel(Hr(k); subplot(222); stem(L,h,filled); axis(0,N-1,-0.1,0.3); xlabel(n ); ylabel(h(n); subplot(223); plot(w/pi,Hr,wL(1:17)/pi,Hrs(1:17); axis(0,1,-0.2,1.2); xlabel(omega (pi); ylabel(Hr(omega); subplot(224); plot(w/pi,20*log10(abs(Hr)/max(abs(Hr); axis(0,1,-50,5);xlabel(omega (pi);ylabel(dB ) grid on; MATLAB 图形为：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 16 页精品资料欢迎下载五、实验心得与体会通过这次实验， 学习了运用窗函数法和频率取样法设计具有线性相位的数字低通、高通、带通、带阻滤波器。同时也从图像上加深了对线性相位的理解和认识。00.5100.51 ( )Hr(k)0102030-0.100.10.20.3nh(n)00.5100.51 ( )Hr()00.51-40-200 ( )dB精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 16 页，共 16 页