用凯泽窗设计线性相位带通FIR滤波

中国地质大学（北京）本科课程报告2018年1月10日中国地质大学（北京）课程名称数字信号处理实验名称用凯泽窗设计线性相位带通FIR滤波学生姓名衡星院（系）地球物理与信息技术学院专业测控技术与仪器学号1010152213指导教师陈玉东目录1实验目的...................................................................................................................-1-2实验环境...................................................................................................................-1-3实验原理及设计流程图...........................................................................................-1-4实验内容及程序源代码...........................................................................................-2-5实验结果...................................................................................................................-5-6实验思考...................................................................................................................-6-7心得与体会...............................................................................................................-7-姓名：衡星学号：1010152213成绩：11实验目的基于MATLAB环境，熟悉窗函数法设计线性相位FIR滤波器的原理和方法。2实验环境硬件环境:计算机,软件环境：MATLAB平台。3实验原理及设计流程图(1)利用凯泽窗原理求解滤波器长度N和形状参数β：凯泽窗是对于给定的阻带衰减，使主瓣具有最大能量意义下的最佳窗函数，因此具有最陡的过渡带。窗函数的形式为：10)(])121(1[)(020NnINnInwk，，其中为形状参数，I0[x]是第一类修正零阶贝塞尔函数，其幂函数展开式为：120])2(!1[1][kkxkxI.若给定滤波器的过渡带宽度(rad)和阻带衰减As(dB)，则滤波器的长度和形状参数可由下列经验公式给出：Δ286.295.7sAN，0.40.1102(8.7)50dB0.5842(21)0.07886(21)2150dB021dBssssssAAAAAA，，，，代入11110.15spsp，As=60dB得N=50，β=5.65.(2)理想带通滤波器单位抽样响应hd(n)的计算：若Hd(ejw)为理想带通线性相位滤波器，12(e)0jccjdeH，，其它，那么，该理想带通滤波器可以根据两个幅值相等而截止频率不同的理想低通滤波器形成，截止频率分别近似为：11122211();()22cpscps，故2121()()21sin[()]sin[()]11()dd22()()ccccjnjnccdnnhneenn，姓名：衡星学号：1010152213成绩：2其中α为线性相位所必须的移位，应满足α=(N-1)/2.(3)设计流程图：计算h(n)=hd(n)·wk(n)用stem(wk)显示wk(n)输入带通指标：p1、p2、s1、s2、As、编写函数dbpfilter_FIR.m用以计算带通滤波器单位抽样响应h(n)计算N(N为奇数)及编写函数kaiser_WF.m计算凯泽窗函数wk(n)计算滤波器的幅频及对数幅频特性并画图打印结果、验证指标要求调用ideal_LP计算hd(n)4实验内容及程序源代码(1)凯泽窗函数文件kaiser_WF.m：functionwk=kaiser_WF(N,beta)%N---凯泽窗序列的长度%beta---凯泽窗参数%wk---凯泽窗序列n=[0:1:N-1];fori=1:1:Ntemp=1-(1-2*n(i)/(N-1))^2;iftemp-eps;disp('负数开方');endbessel_w(i)=bessel_IM(beta*sqrt(temp));endwk=bessel_w/bessel_IM(beta);(2)贝塞尔函数bessel_IM.m：functions=bessel_IM(x)eps=10^(-12);k=1;s=1;D=1;whileD(eps*s)T=x/(2*k);姓名：衡星学号：1010152213成绩：3D=D*T^2;s=s+D;k=k+1;end所用的公式为:120])2(!1[1][kkxkxI.(3)理想低通滤波器单位抽样响应函数ideal_LP.m：function[hd]=ideal_LP(wc,N)%hd---理想低通滤波器单位冲激响应（0=n=(n-1))%wc---截至频率（单位弧度/秒）%N---理想低通滤波器的长度alpha=(N-1)/2;n=[0:(N-1)];m=n-alpha+eps;hd=sin(wc*m)./(pi*m);(4)dbpfilter.m主程序文件：通带截止频率p1=0.3，p2=0.5;s1=0.15，s2=0.65，阻带衰减大于60dB.%d_omega---输入数字频率数组（向量）；%d_omega(1)---阻带下边缘截止频率%d_omega(2)---通带下边缘截止频率%d_omega(3)---通带上边缘截止频率%d_omega(4)---阻带上边缘截止频率%N---数字带通滤波器的长度%hd---理想数字带通滤波器单位冲激响应%h---实际数字带通滤波器单位冲激响应（所设计的）%wk---凯泽窗序列%beta---凯泽窗参数As=60;d_omega=[0.15,0.3,0.5,0.65]*pidelta_w=d_omega(2)-d_omega(1);N=(As-7.95)/(2.286*delta_w)N=ceil(N);M=N;whileM0M=M/2;ifM==ceil(N/2);N=N+1;endendN=N姓名：衡星学号：1010152213成绩：4n=[0:1:N-1];beta=0.1102*(As-8.7)wk=kaiser_WF(N,beta)wc_lower=(d_omega(1)+d_omega(2))/2;wc_upper=(d_omega(3)+d_omega(4))/2;hd=ideal_LP(wc_upper,N)-ideal_LP(wc_lower,N)h=hd.*wk;[H,w]=freqz(h,1,1000,'whole');H=(H(1:501))';w=(w(1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));figure(1);clf;subplot(221);stem(n,hd);title('理想带通滤波器单位冲激响应');axis([-1,N,-0.4,0.6]);ylabel('h_d(n)');subplot(222);stem(n,wk);title('凯泽窗(Kaiser)');axis([-1,N,0,1.3]);ylabel('w_k(n)');subplot(223);stem(n,h);axis([-1,N,-0.4,0.6]);ylabel('h(n)');title('实际带通滤波器单位冲激响应');subplot(224);plot(w/pi,db);axis([0,1,-100,10]);ylabel('分贝');title('实际带通滤波器频率响应(dB)');xlabel('单位为弧度(\pi)');axis([01-8010]);姓名：衡星学号：1010152213成绩：5set(gca,'XTickMode','manual','XTick',d_omega/pi);set(gca,'YTickMode','manual','YTick',[-500]);grid;fileID=fopen('FIR_KAISER_W_h.txt','w');fprintf(fileID,'有限长单位冲激响应滤波器长度N=%d\n\n',N);fprintf(fileID,'%12s\n','h');fprintf(fileID,'%12.8f\n',h);fclose(fileID);5实验结果运行结果如下：d_omega=0.47120.94251.57082.0420N=48.3174N=49beta=5.6533wk=1至13列0.02040.03600.05600.08060.11020.14470.18420.22850.27740.33020.38640.44540.506114至26列0.56770.62900.68910.74670.80070.85010.89380.93090.96060.98230.99561.00000.995627至39列0.98230.96060.93090.89380.85010.80070.74670.68910.62900.56770.50610.44540.386440至49列0.33020.27740.22850.18420.14470.11020.08060.05600.03600.0204wc_lower=0.7069wc_upper=1.8064姓名：衡星学号：1010152213成绩：6hd=1至13列0.0048-0.00180.0106-0.0080-0.0318-0.00880.0130-0.00240.00720.03920.0139-0.0301-0.013314至26列-0.0042-0.0450-0.02130.06120.0478-0.00510.04870.0398-0.1711-0.22950.10280.35000.102827至39列-0.2295-0.17110.03980.0487-0.00510.04780.0612-0.0213-0.0450-0.0042-0.0133-0.03010.013940至49列0.03920.0072-0.00240.0130-0.0088-0.0318-0.00800.0106-0.00180.0048图用凯泽窗设计出的线性相位带通滤波器6实验思考(1)窗函数有哪些指标要求?窗函数有截短和平滑的作用，窗函数选择的好，可以在相同阶次的情况下，提高滤波器的性能，或是在满足设计要求的情况下，减少滤波器阶数。选窗标准：较低的旁瓣幅度，尤其是第一旁瓣；旁瓣幅度要下降得快，以利于增加阻带衰减；主瓣宽度要窄，这样滤波器过渡带较窄。但这三点难以同时满足，当选用主瓣宽度较窄时，虽然得到的幅频特性较陡峭，但通带、阻带波动会明显增加；当选用较低的旁瓣幅度时，虽然得到的幅频特性较平缓匀滑，但过渡带变宽。因此，实际的选择往往取折中。姓名：衡星学号：1010152213成绩：7一般情况下，从“矩形窗、三角窗、汉宁窗、海明窗、布拉克曼窗、凯泽窗”中选择，