利用matlab设计巴特沃斯低通滤波器

三峡大学课程设计报告专业班级 20091421 课程 数字信号处理课程设计学号 学生姓名 姜祥奔 指导教师 2012年5月平时成绩（20%）报告成绩（40%）答辩成绩（40%）总成绩

数字信号处理课程设计实验一:用双线性变换法和脉冲响应不变法设计IIR数字滤波器采用双线性变换法设计一个巴特沃斯数字低通滤波器，要求：通带截止频率100Hz，阻带截止频率200Hz，通带衰减指标Rp小于2dB，阻带衰减Rs大于15dB，滤波器采样频率Fs=500Hz。绘制频率响应曲线。理论部分：(原理及设计过程)第一步：论/和①的关系及数字域性能的公式表示。模拟频率与数字频率之间为1线性关系，①=QT=2时,T=—=10-4,T为抽样周期，故sf=1kHz对应于①=2兀x100x-1-=0.4兀修正后2w、。=tan(^)1T2f=1kHz对应于w=2兀x200x-00=0.8兀修正后O=Ztan(二)1T2按衰减的定义20log10H（e20log10H（ej0）H（ej0.4兀）20log10H（ej0）H（ej0.8兀）>15设w=0处频率响应幅度归一化为1,即H(ej0)=1，则上两式变成20log10H（ej0.4兀）>—2 （1）20log10H（ej0.8兀）<—15 （2）

这就是数字滤波器的性能指标的表达式。2下面把数字低通滤波器的性能要求转变为“样本”模拟低通滤波器的性能要求。由s=OT，按修正式H（ee）wH（jp，设没有混叠效应（即混叠效应设计完成后再进行校验）则有（3）H(ee)=Ha(j?)=七(j。)， ^<（3）利用（3）式，由（1）、（2）式可写出模拟低通滤波器的指标为指标，故应选取比求出的N大一点的整数，故取N=2，代入通带条件(6-118a)式，得Q，c=1.038x103,即20logi0H0.4兀()=20logi0H(j2兀x102)>-2（420logi0H0.4兀()=20logi0H(j2兀x102)>-2（4）20logi0H,.0.8兀、()=20logi0七(j4kx102)<-15...(5)3计算“样本”模拟低通滤波器所需的阶数N及3dB截止频率。C。巴特沃思低通滤波器的幅度平方函数是Ha(j。)|1。1+( )2N。C以分贝形式表示上式，即20log10Hagm-10%。1+( )2N。C把求出的性能指标关系（4）式、（5）式代入（6）式得-10log10，，O、1+（卜）2N。C>-2（6）-10log10。1+(—)2N<-15。C」先用等号来满足指标，可得1+(¥)2N=100.2C=101.5解此两方程，得N=1.3709,N是滤波器阶次，必须取整数，为了满足或超过给定N=2,Q^=1.038x103 (7)显然，利用(7)式的N和巳,阻带指标正好满足，而通带指标则低于2dB,这正好对减少冲激响应不变法造成的频谱混叠效应是有利的。程序部分：wp=100*2*pi;程序部分：wp=100*2*pi;止频率ws=200*2*pi;Rp=2;As=15;标%转换为模拟滤波器指标Fs=500;T=1/Fs;wp1=wp*T;ws1=ws*T;Omgp=(2/T)*tan(wp1/2);Omgs=(2/T)*tan(ws1/2);%模拟原型滤波器计算%数字滤波器的通带截%数字滤波器的阻带截止频率%输入滤波器的通阻带衰减指%原型通带频率预修正%原型通阻频率预修正[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s'); %计算阶数n和截止频率[z0,p0,k0]=buttap(n); %归一化原型设计[ba,aa]=zp2tf(z0,p0,k0); %把滤波器零极点转化为传递函数模型[ba1,aa1]=lp2lp(ba,aa,Omgc); %变换为模拟低通滤波器系数b,a%或者[ba1,aa1]=butter(n,Omgc,’s’);%用双线性变换法计算数字滤波器系数[bd,ad]=bilinear(ba1,aa1,Fs); %用双线性变换法求数字滤波器系数b,a%求数字系统的频率特性[H,w]=freqz(bd,ad);plot(w*Fs/(2*pi),20*log10(abs(H)));ylabel('增益/dB');xlabel('频率(Hz)');title('数字滤波器幅度响应');gridon;数字滤波器的频率响应图:数字滤波器幅度响应-100 1 1 1 L_□ 数字滤波器幅度响应-100 1 1 1 L_□ 50 100 150 200OOOOOOOOO123456789---------250实验二：用窗函数法设计FIR^字^波器选择合适的窗函数设计一个FIR数字低通滤波器，要求：带通截止频率为Wp=0.2n,阻带截止频率Ws=0.3n,阻带衰减不小于40dB,通带衰减不大于3dB。描绘滤波器的幅频响应曲线。理论部分：(原理及设计过程)1.数字频率。通带截止频率:w通带截止频率:w=0.2兀

p阻带截止频率:阻带衰减：32=4迦2.求气(n)。设Hd(ejw)为理想线性相位滤波器H(ejw)首先由所需低通滤波器的过渡带求理想低通滤波器的截止频率Q(由于Q为两C C

个肩峰值处的频率中点，而由Q〃到Q,之间的过渡带宽并非两个肩峰间的频率差，因而以下求出的Qc有一定的近似)，并将其转化为对应的数字频率为w=2(w+w,)=0.25兀由此可得h(n)=—!—jwh(n)=—!—jwcejw(n-t)dw=<d2兀-wc sinlwc兀(n-t)w—c，冗(n-t)]n方n=兀其中T为线性相位所必需的移位，且T求窗函数。由阻带衰减52=40dB，查课本中的表7-3可选汉宁窗，其阻带最小衰减44dB满足要求。所要求的过度带宽(数字频域)Aw=w-w=0.1兀p由于汉宁窗过渡带宽满足侦=6N所以TOC\o"1-5"\h\z6.2兀 〜N= =62AwN-1平t= =312求h(n)。由汉宁窗表达式w(n)确定FIR滤波器的h(n)。汉宁窗w(n)=2[1-cos(Nn^R^n.N-1sin[w(n-—2—)]汉宁窗hd(n)= N-1兀(n )2所以h(n)所以h(n)=h^(n)*w(n)由h(n)求得H(e]W)，并检验各项指标是否满足要求。如不满足要求，则要改变N,或改变窗形状(或两者都改变)，然后重新计算。\

程序部分:子函数1:functionhd=ideal_lp(wc,M);alpha=(M-1)/2;n=[0:1:(M-1)];m=n-alpha+eps; %addsmallestnumbertoavoidividedbyzerohd=sin(wc*m)./(pi*m);子函数2：function[db,mag,pha,w]=freqz_m4(b,a)[H,w]=freqz(b,a,1000,'whole');function[db,mag,pha,w]=freqz_m4(b,a)[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501));w=(w(1:1:501));mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);主程序：clearallWp=0.2*pi;Ws=0.3*pi;tr_width=Ws-Wp;N=ceil(6.2*pi/tr_width)+1;n=0:1:N-1;Wc=(Ws+Wp)/2;hd=ideal_lp(Wc,N);w_ham=(hann(N))';h=hd.*w_ham;[db,mag,pha,w]=freqz_m4(h,[1]);delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1)));As=-round(max(db(Ws/delta_w+1:1:501)));subplot(221)stem(n,hd)title('理想单位脉冲响应hd(n)')subplot(222)stem(n,w_ham)title('汉宁窗w(n)')subplot(223)stem(n,h)title('实际单位脉冲响应h(n)')%过渡带宽度%滤波器长度%理想低通滤波器的截止频率%理想低通滤波器的单位冲激响应%汉宁窗%截取得到实际的单位脉冲响应%计算实际滤波器的幅度响应%实际通带纹波%实际阻带纹波%火柴棍图subplot(224)plot(w/pi,db)title('幅度响应(dB)')axis([0,1,-100,10])FIR低通滤波器特性图:不得不承认，这一次数字信号处理的课程设计过程是艰苦的，但同时又是收获丰富的。因为我们是上学期修的数字信号处理，而课程设计却安排在下学期，经过半年的时间，很多学过的东西都变得模糊的。刚选定课题时，自己甚至连题目的意思都看不懂，“什么是双线性变换法？”“什么是巴特沃思数字低通滤波器？”，一系列的问号在头脑中打转。Matlab程序设计经验的缺乏也曾一度令自己失去持续做程题的勇气。现在，虽然不能拍着胸脯说