实验四窗函数法设计FIR数字滤波器.doc

word实验四 窗函数法设计FIR数字滤波器一、实验目的1、掌握窗函数法设计FIR数字滤波器的原理与具体方法。 2、掌握频率取样法设计FIR数字滤波器的原理和根本方法。 3、学习利用窗函数法和频率取样法设计低通、带通、高通、带阻数字滤波器。二、实验环境计算机、MATLAB软件三、实验根底理论窗函数设计FIR滤波器窗函数设计法的根本思想为，首先选择一个适当的理想的滤波器，然后用窗函数截取它的单位脉冲响应，得到线性相位和因果的FIR滤波器。这种方法的重点是选择一个适宜的窗函数和理想滤波器，使设计的滤波器的单位脉冲响应逼近理想滤波器的单位脉冲响应。1给定理想滤波器的频率响应，在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽为的低通滤波器由下式给定：其中为采样延迟，其作用是为了得到一个因果系统。2确定这个滤波器的单位脉冲响应为了得到一个长度为N的因果的线性相位FIR滤波器，我们令3用窗函数截取得到所设计FIR数字滤波器：常用的窗函数有矩形Rectangular窗，汉宁Hanning窗，海明Hamming窗、布莱克曼Blackman窗、凯瑟Kaiser窗等表4-1 MATLAB中产生窗函数的命令MATLAB函数窗函数MATLAB函数窗函数Boxcar矩形窗函数Blackman布莱克曼窗Hanning汉宁窗函数Kaiser凯瑟窗函数Hamming海明窗表4-2 常用窗函数的特性窗函数窗函数频率特性加窗后滤波器指标旁瓣峰值dB主瓣宽度过渡带宽最小阻带衰减dB矩形窗-134/N/N-21汉宁窗-318/N/N-44海明窗-418/N/N-53布莱克曼窗-5712/N11/N-74凯瑟窗是一种广泛在实际中广泛应用的窗函数，它由下式给定：其中是修正的零阶贝塞尔函数，参数控制最小阻带衰减，这种窗函数对于一样的N可以提供不同的过渡带宽。由于贝塞尔函数比拟复杂，这种窗函数的设计方程很难推导，然而幸运的是，有一些经验设计方程可以直接使用。给定的指标，滤波器长度N和凯瑟窗参数可以按如下凯瑟窗方程给出过渡带带宽：频率取样设计FIR滤波器频率取样法从频域出发，把理想的滤波器等间隔采样得到,将作为实际设计滤波器的：得到以后可以由来确定唯一确定滤波器的单位脉冲响应,可以由求得：其中为内插函数：有求得的频率响应将逼近。如果我们设计的是线性相位FIR滤波器，如此的幅度和相位满足线性相位滤波器的约束条件。我们将表示为如下形式当为实数，如此由此得到即为中心偶对称。在利用线性相位条件可知，对于1型和2型线性相位滤波器：对于3型和4型线性相位滤波器1由给定的理想滤波器给出和。2由求得3根据求得或四、实验内容1、设计一个数字低通FIR滤波器，其技术指标如下：分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。（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图形：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);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');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)');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');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滤波器，其技术指标如下：下阻带边缘：下通带边缘：上通带边缘：上阻带边缘：程序代码：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;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数字低通滤波器，满足以下指标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图形如下：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');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);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图形:各阶数的通阻带指标如下表：阶数阻带衰减50dB是否满足指标N=20不满足不满足否N=40不满足满足否N=60满足满足是对于高通滤波器，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图形为：五、实验心得与体会通过这次实验，学习了运用窗函数法和频率取样法设计具有线性相位的数字低通、高通、带通、带阻滤波器。同时也从图像上加深了对线性相位的理解和认识。文案大全